CMIP6 prep

---

library(raster)

library(plyr)
library(dplyr)

library(ggplot2)
library(viridisLite)
library(colorspace)
library(RColorBrewer)

library(sf)
library(stringr)

library(ggpubr)
library(kableExtra)

memory.limit(size = 160000)
## [1] Inf

#display.brewer.pal(n = 11, name ="RdYlGn")
colorvec <- brewer.pal(n = 11, name ="RdYlGn")

1 Data

land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land


elev <- raster("Worldclim/wc2.1_10m/wc2.1_10m_elev.tif") 
elev.df <- elev %>% projectRaster(to = land_mask) %>% as.data.frame(xy = T) %>% setNames(c("lon","lat", "z"))


ipcc_regions <- shapefile("Masks/IPCC-WGI-reference-regions-v4.shp") %>% spTransform(crs("EPSG:4326")) 
ipcc_regions.raster <- ipcc_regions %>% rasterize(land_mask)
ipcc_regions.df <- as.data.frame(ipcc_regions.raster, xy = T) %>% setNames(c("lon","lat","Continent","Type","Name","Acronym"))

tas.all_annual <- read.table("CMIP6/tas/tas.all_annual.txt")
pr.all_annual <- read.table("CMIP6/pr/pr.all_annual.txt")
sfcWind.all_annual <- read.table("CMIP6/sfcWind/sfcWind.all_annual.txt")
hfls.all_annual <- read.table("CMIP6/hfls/hfls.all_annual.txt")
hfss.all_annual <- read.table("CMIP6/hfss/hfss.all_annual.txt")
hurs.all_annual <- read.table("CMIP6/hurs/hurs.all_annual.txt")

cmip6_annual <- tas.all_annual %>% merge(pr.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_annual, "cmip6_annual.txt")

tas.all_january <- read.table("CMIP6/tas/tas.all_january.txt")
pr.all_january <- read.table("CMIP6/pr/pr.all_january.txt")
sfcWind.all_january <- read.table("CMIP6/sfcWind/sfcWind.all_january.txt")
hfls.all_january <- read.table("CMIP6/hfls/hfls.all_january.txt")
hfss.all_january <- read.table("CMIP6/hfss/hfss.all_january.txt")
hurs.all_january <- read.table("CMIP6/hurs/hurs.all_january.txt")

cmip6_january <- tas.all_january %>% merge(pr.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_january, "cmip6_january.txt")

tas.all_february <- read.table("CMIP6/tas/tas.all_february.txt")
pr.all_february <- read.table("CMIP6/pr/pr.all_february.txt")
sfcWind.all_february <- read.table("CMIP6/sfcWind/sfcWind.all_february.txt")
hfls.all_february <- read.table("CMIP6/hfls/hfls.all_february.txt")
hfss.all_february <- read.table("CMIP6/hfss/hfss.all_february.txt")
hurs.all_february <- read.table("CMIP6/hurs/hurs.all_february.txt")

cmip6_february <- tas.all_february %>% merge(pr.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_february, "cmip6_february.txt")

tas.all_march <- read.table("CMIP6/tas/tas.all_march.txt")
pr.all_march <- read.table("CMIP6/pr/pr.all_march.txt")
sfcWind.all_march <- read.table("CMIP6/sfcWind/sfcWind.all_march.txt")
hfls.all_march <- read.table("CMIP6/hfls/hfls.all_march.txt")
hfss.all_march <- read.table("CMIP6/hfss/hfss.all_march.txt")
hurs.all_march <- read.table("CMIP6/hurs/hurs.all_march.txt")

cmip6_march <- tas.all_march %>% merge(pr.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_march, "cmip6_march.txt")

tas.all_april <- read.table("CMIP6/tas/tas.all_april.txt")
pr.all_april <- read.table("CMIP6/pr/pr.all_april.txt")
sfcWind.all_april <- read.table("CMIP6/sfcWind/sfcWind.all_april.txt")
hfls.all_april <- read.table("CMIP6/hfls/hfls.all_april.txt")
hfss.all_april <- read.table("CMIP6/hfss/hfss.all_april.txt")
hurs.all_april <- read.table("CMIP6/hurs/hurs.all_april.txt")

cmip6_april <- tas.all_april %>% merge(pr.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_april, "cmip6_april.txt")

tas.all_may <- read.table("CMIP6/tas/tas.all_may.txt")
pr.all_may <- read.table("CMIP6/pr/pr.all_may.txt")
sfcWind.all_may <- read.table("CMIP6/sfcWind/sfcWind.all_may.txt")
hfls.all_may <- read.table("CMIP6/hfls/hfls.all_may.txt")
hfss.all_may <- read.table("CMIP6/hfss/hfss.all_may.txt")
hurs.all_may <- read.table("CMIP6/hurs/hurs.all_may.txt")

cmip6_may <- tas.all_may %>% merge(pr.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_may, "cmip6_may.txt")

tas.all_june <- read.table("CMIP6/tas/tas.all_june.txt")
pr.all_june <- read.table("CMIP6/pr/pr.all_june.txt")
sfcWind.all_june <- read.table("CMIP6/sfcWind/sfcWind.all_june.txt")
hfls.all_june <- read.table("CMIP6/hfls/hfls.all_june.txt")
hfss.all_june <- read.table("CMIP6/hfss/hfss.all_june.txt")
hurs.all_june <- read.table("CMIP6/hurs/hurs.all_june.txt")

cmip6_june <- tas.all_june %>% merge(pr.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_june, "cmip6_june.txt")

tas.all_july <- read.table("CMIP6/tas/tas.all_july.txt")
pr.all_july <- read.table("CMIP6/pr/pr.all_july.txt")
sfcWind.all_july <- read.table("CMIP6/sfcWind/sfcWind.all_july.txt")
hfls.all_july <- read.table("CMIP6/hfls/hfls.all_july.txt")
hfss.all_july <- read.table("CMIP6/hfss/hfss.all_july.txt")
hurs.all_july <- read.table("CMIP6/hurs/hurs.all_july.txt")

cmip6_july <- tas.all_july %>% merge(pr.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_july, "cmip6_july.txt")

tas.all_august <- read.table("CMIP6/tas/tas.all_august.txt")
pr.all_august <- read.table("CMIP6/pr/pr.all_august.txt")
sfcWind.all_august <- read.table("CMIP6/sfcWind/sfcWind.all_august.txt")
hfls.all_august <- read.table("CMIP6/hfls/hfls.all_august.txt")
hfss.all_august <- read.table("CMIP6/hfss/hfss.all_august.txt")
hurs.all_august <- read.table("CMIP6/hurs/hurs.all_august.txt")

cmip6_august <- tas.all_august %>% merge(pr.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_august, "cmip6_august.txt")

tas.all_september <- read.table("CMIP6/tas/tas.all_september.txt")
pr.all_september <- read.table("CMIP6/pr/pr.all_september.txt")
sfcWind.all_september <- read.table("CMIP6/sfcWind/sfcWind.all_september.txt")
hfls.all_september <- read.table("CMIP6/hfls/hfls.all_september.txt")
hfss.all_september <- read.table("CMIP6/hfss/hfss.all_september.txt")
hurs.all_september <- read.table("CMIP6/hurs/hurs.all_september.txt")

cmip6_september <- tas.all_september %>% merge(pr.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_september, "cmip6_september.txt")

tas.all_october <- read.table("CMIP6/tas/tas.all_october.txt")
pr.all_october <- read.table("CMIP6/pr/pr.all_october.txt")
sfcWind.all_october <- read.table("CMIP6/sfcWind/sfcWind.all_october.txt")
hfls.all_october <- read.table("CMIP6/hfls/hfls.all_october.txt")
hfss.all_october <- read.table("CMIP6/hfss/hfss.all_october.txt")
hurs.all_october <- read.table("CMIP6/hurs/hurs.all_october.txt")

cmip6_october <- tas.all_october %>% merge(pr.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_october, "cmip6_october.txt")

tas.all_november <- read.table("CMIP6/tas/tas.all_november.txt")
pr.all_november <- read.table("CMIP6/pr/pr.all_november.txt")
sfcWind.all_november <- read.table("CMIP6/sfcWind/sfcWind.all_november.txt")
hfls.all_november <- read.table("CMIP6/hfls/hfls.all_november.txt")
hfss.all_november <- read.table("CMIP6/hfss/hfss.all_november.txt")
hurs.all_november <- read.table("CMIP6/hurs/hurs.all_november.txt")

cmip6_november <- tas.all_november %>% merge(pr.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_november, "cmip6_november.txt")

tas.all_december <- read.table("CMIP6/tas/tas.all_december.txt")
pr.all_december <- read.table("CMIP6/pr/pr.all_december.txt")
sfcWind.all_december <- read.table("CMIP6/sfcWind/sfcWind.all_december.txt")
hfls.all_december <- read.table("CMIP6/hfls/hfls.all_december.txt")
hfss.all_december <- read.table("CMIP6/hfss/hfss.all_december.txt")
hurs.all_december <- read.table("CMIP6/hurs/hurs.all_december.txt")

cmip6_december <- tas.all_december %>% merge(pr.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_december, "cmip6_december.txt")

2 Aridity Index

2.1 De Martonne

cmip6 <- read.table("cmip6_annual.txt")

breaks.martonne <- c(0,5,10,20,"30","40",Inf)
cat.martonne <- c("[0,5]" = "Desert", "(5,10]"= "Arid", "(10,20]" = "Semi-arid","(20,30]" = "Temperate", "(30,40]" = "Humid", "(40,Inf]" = "Forest")
col.martonne <- c("Desert" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Temperate" = colorvec[8], "Humid"= colorvec[10],"Forest"= colorvec[11])

cmip6$AIm <- with(cmip6,pr*60*60*24*365/(tas-273.15+10)) # pr in mm/y, tas in Ceslsius


cmip6$cat.AIm <- cmip6$AIm %>% cut(breaks.martonne, include.lowest = T) %>% revalue(cat.martonne)


map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1970_2000" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AIm))+
    scale_fill_manual(values = col.martonne, na.translate = F)+
    labs(title = paste(i, "1970-2000", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

2.2 UNEP

ea from Rh:

ea = 6.108(RH/100)exp(17.27*T/(T+237.3)) Refs https://www.nature.com/articles/s41598-023-40499-6 : Allen and Tetens

id_vars <- c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")

cmip6_annual <- mutate(read.table("cmip6_annual.txt"),
                       t = tas - 273.15,
                       Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                       P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                       es = 0.61078*exp(17.27*t/(t+237.3)),
                       ea = (hurs/100)*es,
                       delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                       gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                       ET0_annual = ((0.408*delta*Rn+gamma*900/((t)+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind)))) 

cmip6_01 <- mutate(read.table("cmip6_january.txt"),
                   t = tas - 273.15,
                   pr_01 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_01 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind)))) # Evapotranspiration, G considered to be negligible. in mm/day. Wind speed at 10 m is converted to wind speed at 2 m by multiplying by 0.748

cmip6_02 <- mutate(read.table("cmip6_february.txt"),
                   t = tas - 273.15,
                   pr_02 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_02 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_03 <- mutate(read.table("cmip6_march.txt"),
                   t = tas - 273.15,
                   pr_03 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_03 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_04 <- mutate(read.table("cmip6_april.txt"),
                   t = tas - 273.15,
                   pr_04 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_04 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_05 <- mutate(read.table("cmip6_may.txt"),
                   t = tas - 273.15,
                   pr_05 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_05 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_06 <- mutate(read.table("cmip6_june.txt"),
                   t = tas - 273.15,
                   pr_06 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_06 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_07 <- mutate(read.table("cmip6_july.txt"),
                   t = tas - 273.15,
                   pr_07 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_07 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_08 <- mutate(read.table("cmip6_august.txt"),
                   t = tas - 273.15,
                   pr_08 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_08 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_09 <- mutate(read.table("cmip6_september.txt"),
                   t = tas - 273.15,
                   pr_09 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_09 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_10 <- mutate(read.table("cmip6_october.txt"),
                   t = tas - 273.15,
                   pr_10 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_10 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_11 <- mutate(read.table("cmip6_november.txt"),
                   t = tas - 273.15,
                   pr_11 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_11 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_12 <- mutate(read.table("cmip6_december.txt"),
                   t = tas - 273.15,
                   pr_12 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_12 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

et0_monthly <- select(cmip6_annual, id_vars) %>%
  merge(select(cmip6_01, c(id_vars, "ET0_01", "pr_01")), by = id_vars) %>%
  merge(select(cmip6_02, c(id_vars, "ET0_02", "pr_02")), by = id_vars) %>%
  merge(select(cmip6_03, c(id_vars, "ET0_03", "pr_03")), by = id_vars) %>%
  merge(select(cmip6_04, c(id_vars, "ET0_04", "pr_04")), by = id_vars) %>%
  merge(select(cmip6_05, c(id_vars, "ET0_05", "pr_05")), by = id_vars) %>%
  merge(select(cmip6_06, c(id_vars, "ET0_06", "pr_06")), by = id_vars) %>%
  merge(select(cmip6_07, c(id_vars, "ET0_07", "pr_07")), by = id_vars) %>%
  merge(select(cmip6_08, c(id_vars, "ET0_08", "pr_08")), by = id_vars) %>%
  merge(select(cmip6_09, c(id_vars, "ET0_09", "pr_09")), by = id_vars) %>%
  merge(select(cmip6_10, c(id_vars, "ET0_10", "pr_10")), by = id_vars) %>%   
  merge(select(cmip6_11, c(id_vars, "ET0_11", "pr_11")), by = id_vars) %>%
  merge(select(cmip6_12, c(id_vars, "ET0_12", "pr_12")), by = id_vars)

et0_monthly$spr <- with(et0_monthly, pr_01*31 + pr_02 * 28 + pr_03 * 31 + pr_04 * 30 + pr_05 * 31 + pr_06 * 30 +
                          pr_07 * 31 + pr_08 * 31 + pr_09 * 30 + pr_10 * 31 + pr_11 * 30 + pr_12 * 31)
et0_monthly$sET0 <- with(et0_monthly, ET0_01*31 + ET0_02 * 28 + ET0_03 * 31 + ET0_04 * 30 + ET0_05 * 31 + ET0_06 * 30 +
                           ET0_07 * 31 + ET0_08 * 31 + ET0_09 * 30 + ET0_10 * 31 + ET0_11 * 30 + ET0_12 * 31)

cmip6 <- merge(cmip6_annual, et0_monthly, by = id_vars)
cmip6$AI <- with(cmip6, abs(spr/sET0))  
cmip6$AI_annual <- with(cmip6, abs((pr*60*60*24*365)/(ET0_annual*365)))

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6$cat.AI[which(cmip6$sET0 < 400)] <- "Cold"
cmip6$cat.AI %>% unique()

cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>% 
  revalue(cat.unesco)


cmip6$cat.AI_annual <- cmip6$AI_annual %>% cut(breaks = breaks.unesco) %>% 
  revalue(cat.unesco)

cmip6$cat.AI_annual[which(cmip6$ET0_annual*365 < 400)] <- "Cold"

write.table(cmip6, "cmip6_AI.txt")

cmip6 <- read.table("cmip6_AI.txt")

cmip6_AI <- cmip6 %>% subset(period == "1970_2000") %>%
  group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, z) %>% 
  dplyr::summarise(tas = mean(tas, na.rm = T) - 273.15, pr = mean(pr, na.rm = T) * 60 * 60 * 24 * 365,
                   sfcWind = mean(sfcWind, na.rm = T), Rn = mean(Rn, na.rm = T),
                   ea = mean(ea, na.rm = T), ET0 = mean(sET0, na.rm = T), AI = mean(AI, na.rm = T)) %>%
  ungroup() %>%  
  mutate(source = "CMIP6")%>%
  select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI"))


write.table(cmip6_AI, "cmip6ref.txt")

3 Maps AI by model

cmip6 <- read.table("cmip6_AI.txt")

breaks.unesco <- c(-Inf,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(-Inf,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid")

col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

3.1 1850-1880

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1850_1880" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1850-1880", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "1850_1880" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))


tab <- cmip6 %>% subset(period == "1850_1880" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	3.6	6.3	9.2	4.3	23.4	26.4	47.7	2.6
CESM	3.6	7.5	10.2	3.5	24.8	24.8	47.8	2.6
CMCC	3.9	5.0	6.6	3.7	19.2	29.3	48.9	2.6
CMCC-ESM2	4.2	5.7	5.2	1.9	17.0	16.7	63.6	2.6
CNRM	4.0	8.4	8.0	4.6	25.0	26.4	46.0	2.6
EC-Earth3	6.0	6.2	7.6	3.1	22.9	25.6	48.9	2.6
FGOALS	3.9	8.9	11.0	5.0	28.8	20.8	47.8	2.6
GFDL-ESM4	3.9	6.5	6.1	2.9	19.4	25.8	52.3	2.6
INM	3.0	5.7	9.4	4.0	22.1	28.2	47.0	2.6
INM-CM5	2.4	6.2	8.0	3.9	20.5	30.0	46.9	2.6
MPI	5.9	7.6	7.0	3.0	23.5	24.0	49.9	2.6
MRI	4.4	7.4	6.0	2.7	20.5	27.5	49.5	2.6
NorESM-2-MM	3.6	7.3	12.8	4.1	27.8	24.7	45.0	2.6
mean	4.0	6.8	8.2	3.6	22.6	25.4	49.3	2.6
sd	1.0	1.1	2.2	0.9	5.2	3.5	4.7	0.0

3.2 1970-2000

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1970_2000" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1850-1880", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

library(ggrepel) 

pie <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","POLAR","ATLANTIC","INDIAN","ARCTIC")) %>%
  select(c("source", "cat.AI")) %>%
  reshape2::melt(id.vars = c("source")) %>%
  group_by(source, value) %>%
  summarise(count = n()) %>%
  mutate(perc = round(count/sum(count)*100, 1),
         lab.y = rev(cumsum(rev(count))) - count*0.5) %>% 
  ungroup() 

pie_list <- list()

for(i in unique(pie$source)){
  p <- ggplot(subset(pie, source == i))+
    geom_bar(aes(x = "", y = count, fill = value), width = 1, stat = "identity")+
    coord_polar("y", start = 0)+
    scale_fill_manual(values = col.cat)+
    geom_label_repel(aes(x = "", y = lab.y, label = perc), nudge_x = 0.5)+
    labs(fill = "", title = i)+
    theme_void()+theme(legend.position = "right")
  pie_list[[i]] <- p
}

ggpubr::ggarrange(plotlist = pie_list, nrow = 4, ncol = 4, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))


tab <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	3.6	6.0	9.2	4.4	23.2	26.0	48.2	2.6
CESM	3.7	7.4	9.8	3.6	24.5	24.3	48.5	2.6
CMCC	3.6	5.2	7.0	3.8	19.6	29.4	48.5	2.6
CMCC-ESM2	4.5	5.5	5.3	1.8	17.1	17.0	63.3	2.6
CNRM	4.0	8.5	7.9	4.5	24.9	26.6	45.9	2.6
EC-Earth3	5.9	6.6	7.5	3.1	23.1	24.8	49.6	2.6
FGOALS	4.1	7.6	7.0	3.4	22.1	17.5	57.9	2.6
GFDL-ESM4	4.0	6.2	6.2	2.8	19.2	24.6	53.6	2.6
INM	3.1	5.9	9.2	4.0	22.2	27.9	47.4	2.6
INM-CM5	2.3	5.9	8.4	3.9	20.5	29.9	47.0	2.6
MPI	6.0	7.7	7.1	3.0	23.8	23.7	49.9	2.6
MRI	4.3	7.5	6.4	2.7	20.9	27.0	49.5	2.6
NorESM-2-MM	3.6	7.0	11.9	3.9	26.4	23.3	47.7	2.6
mean	4.1	6.7	7.9	3.5	22.2	24.8	50.5	2.6
sd	1.0	1.0	1.8	0.8	4.6	3.9	5.0	0.0

3.3 1985-2010

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1985_2015" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1985-2010", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "1985_2015" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))


tab <- cmip6 %>% subset(period == "1985_2015" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	3.7	5.9	9.6	4.6	23.8	25.9	47.7	2.6
CESM	3.3	7.5	10.3	3.7	24.8	24.6	47.8	2.6
CMCC	3.6	5.3	7.2	4.0	20.1	29.8	47.4	2.6
CMCC-ESM2	4.3	5.9	5.4	1.7	17.3	17.1	63.0	2.6
CNRM	4.0	8.4	7.9	4.3	24.6	27.3	45.5	2.6
EC-Earth3	5.8	6.8	7.4	3.1	23.1	25.6	48.8	2.6
FGOALS	4.0	7.9	7.2	3.4	22.5	17.6	57.3	2.6
GFDL-ESM4	4.1	6.3	6.3	2.7	19.4	25.4	52.6	2.6
INM	2.8	5.4	7.9	3.6	19.7	25.4	52.3	2.6
INM-CM5	2.4	5.8	8.3	3.8	20.3	28.3	48.8	2.6
MPI	6.1	7.8	7.2	3.1	24.2	23.7	49.6	2.6
MRI	4.3	7.5	6.3	2.6	20.7	26.7	50.0	2.6
NorESM-2-MM	3.3	6.9	11.1	3.9	25.2	21.1	51.1	2.6
mean	4.0	6.7	7.9	3.4	22.0	24.5	50.9	2.6
sd	1.0	1.0	1.6	0.8	4.4	3.8	4.7	0.0

3.4 SSP 2-4.5

3.4.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP245" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP245" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	3.5	6.4	10.2	4.6	24.7	26.5	46.2	2.6
CESM	3.3	8.6	10.6	3.9	26.4	24.8	46.2	2.6
CMCC	3.1	5.9	7.6	4.2	20.8	30.8	45.8	2.6
CMCC-ESM2	4.2	6.3	5.5	2.0	18.0	17.6	61.7	2.6
CNRM	4.2	8.8	8.4	4.7	26.1	27.7	43.6	2.6
EC-Earth3	5.2	6.4	7.7	3.1	22.4	25.3	49.7	2.6
FGOALS	4.8	9.3	12.9	5.4	32.4	21.1	44.0	2.6
GFDL-ESM4	3.8	7.1	6.6	3.1	20.6	25.7	51.3	2.6
INM	3.2	5.8	9.6	4.1	22.7	28.1	46.5	2.6
INM-CM5	2.5	6.0	8.5	3.8	20.8	30.0	46.7	2.6
MPI	6.2	8.3	7.7	3.1	25.3	24.5	47.6	2.6
MRI	4.3	7.8	6.5	2.8	21.4	27.8	48.1	2.6
NorESM-2-MM	3.7	8.2	11.7	3.9	27.5	23.3	46.5	2.6
mean	4.0	7.3	8.7	3.7	23.7	25.6	48.0	2.6
sd	1.0	1.2	2.2	0.9	5.3	3.6	4.6	0.0

3.4.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP245" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP245" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	3.3	7.3	10.4	5.1	26.1	26.6	44.8	2.6
CESM	4.2	8.1	10.7	3.7	26.7	25.7	44.9	2.6
CMCC	3.0	5.6	8.5	4.5	21.6	32.4	43.4	2.6
CMCC-ESM2	4.3	6.6	5.9	1.9	18.7	18.7	60.0	2.6
CNRM	4.4	9.0	8.6	4.7	26.7	28.5	42.2	2.6
EC-Earth3	5.3	6.4	7.7	3.0	22.4	26.0	49.1	2.6
FGOALS	4.8	9.8	12.6	5.6	32.8	21.1	43.5	2.6
GFDL-ESM4	4.1	6.7	6.8	3.1	20.7	26.8	50.0	2.6
INM	3.3	5.7	9.8	4.1	22.9	28.3	46.3	2.6
INM-CM5	2.6	6.6	7.8	4.1	21.1	30.3	46.0	2.6
MPI	6.2	8.3	7.7	3.1	25.3	24.5	47.6	2.6
MRI	4.5	7.6	6.7	2.8	21.6	28.4	47.4	2.6
NorESM-2-MM	4.1	8.2	11.8	4.1	28.2	23.7	45.6	2.6
mean	4.2	7.4	8.8	3.8	24.2	26.2	47.0	2.6
sd	1.0	1.3	2.1	1.0	5.4	3.7	4.5	0.0

3.5 SSP 3-7.0

3.5.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP370" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP370" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	3.4	6.7	9.9	4.8	24.8	26.0	46.6	2.6
CESM	3.6	7.9	10.8	3.7	26.0	24.8	46.4	2.6
CMCC	2.9	5.8	7.9	4.1	20.7	30.7	46.0	2.6
CMCC-ESM2	2.3	5.8	6.2	2.5	16.8	18.7	61.8	2.6
CNRM	4.3	9.0	8.9	5.0	27.2	27.6	42.6	2.6
EC-Earth3	5.5	6.8	8.0	3.3	23.6	26.8	47.0	2.6
FGOALS	4.2	8.3	8.3	3.7	24.5	19.2	53.8	2.6
GFDL-ESM4	4.4	6.6	6.9	2.9	20.8	25.0	51.8	2.6
INM	3.2	5.8	9.4	4.1	22.5	28.0	46.9	2.6
INM-CM5	2.4	6.2	8.4	3.9	20.9	30.1	46.4	2.6
MPI	6.1	8.5	7.6	3.1	25.3	25.3	46.8	2.6
MRI	4.5	7.5	6.5	2.5	21.0	28.0	48.3	2.6
NorESM-2-MM	3.8	7.6	12.5	4.0	27.9	23.1	46.4	2.6
mean	3.9	7.1	8.6	3.7	23.3	25.6	48.5	2.6
sd	1.1	1.1	1.8	0.8	4.8	3.6	4.8	0.0

3.5.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP370" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP370" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	3.8	7.3	10.7	4.8	26.6	26.0	44.9	2.6
CESM	4.2	8.5	10.8	4.0	27.5	26.0	43.9	2.6
CMCC	3.6	6.3	8.7	5.1	23.7	30.7	43.0	2.6
CMCC-ESM2	4.8	6.9	5.7	2.1	19.5	18.0	59.9	2.6
CNRM	5.1	9.0	9.2	4.7	28.0	27.3	42.2	2.6
EC-Earth3	5.3	6.9	8.6	3.3	24.1	28.1	45.3	2.6
FGOALS	4.7	9.3	9.2	3.0	26.2	19.5	51.8	2.6
GFDL-ESM4	5.2	5.9	7.6	3.3	22.0	25.7	49.7	2.6
INM	4.2	5.7	10.2	4.3	24.4	27.3	45.7	2.6
INM-CM5	2.7	7.2	9.8	4.3	24.0	28.1	45.3	2.6
MPI	6.1	8.5	7.6	3.1	25.3	25.3	46.8	2.6
MRI	5.0	8.3	6.9	2.9	23.1	27.2	47.2	2.6
NorESM-2-MM	4.1	8.2	12.2	4.1	28.6	23.5	45.3	2.6
mean	4.5	7.5	9.0	3.8	24.8	25.6	47.0	2.6
sd	0.9	1.2	1.8	0.9	4.8	3.5	4.7	0.0

3.6 SSP 5-8.5

3.6.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP585" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP585" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	3.4	6.4	10.5	5.2	25.5	26.2	45.8	2.6
CESM	4.0	8.0	10.9	3.6	26.5	25.0	45.9	2.6
CMCC	3.0	5.6	8.0	4.5	21.1	31.2	45.0	2.6
CMCC-ESM2	4.5	6.5	5.6	1.7	18.3	17.7	61.4	2.6
CNRM	4.3	9.2	8.6	4.8	26.9	27.8	42.6	2.6
EC-Earth3	5.4	6.6	8.2	3.2	23.4	27.0	47.0	2.6
FGOALS	4.8	9.4	12.6	5.4	32.2	21.4	43.9	2.6
GFDL-ESM4	4.3	6.6	6.8	3.3	21.0	25.3	51.2	2.6
INM	3.2	5.8	9.7	4.1	22.8	28.3	46.2	2.6
INM-CM5	2.5	6.0	8.6	4.0	21.1	30.3	46.0	2.6
MPI	6.0	8.7	7.7	3.1	25.5	25.7	46.2	2.6
MRI	4.2	7.7	6.7	2.8	21.4	28.4	47.6	2.6
NorESM-2-MM	3.8	7.9	12.2	4.1	28.0	23.1	46.2	2.6
mean	4.1	7.3	8.9	3.8	24.1	26.0	47.3	2.6
sd	1.0	1.3	2.1	1.0	5.4	3.7	4.7	0.0

3.6.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP585" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP585" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(13,14), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	4.3	7.5	11.1	4.5	27.4	26.1	43.9	2.6
CMCC	3.1	6.5	9.4	4.6	23.6	32.5	41.3	2.6
CMCC-ESM2	3.9	7.3	6.2	2.0	19.4	19.5	58.5	2.6
CNRM	4.2	9.4	8.9	4.4	26.9	28.5	42.0	2.6
EC-Earth3	5.4	6.7	9.1	3.4	24.6	28.2	44.6	2.6
FGOALS	10.8	10.4	12.7	4.7	38.6	16.9	41.8	2.6
GFDL-ESM4	4.5	6.9	7.5	3.3	22.2	26.2	49.0	2.6
INM	3.1	6.1	10.3	4.1	23.6	29.0	44.8	2.6
INM-CM5	2.7	6.8	8.3	4.2	22.0	30.7	44.6	2.6
MPI	6.1	8.6	7.7	3.3	25.7	25.6	46.2	2.6
MRI	4.9	8.0	7.4	2.8	23.1	27.9	46.4	2.6
NorESM-2-MM	4.4	8.9	12.1	4.5	29.9	23.9	43.5	2.6
mean	4.8	7.8	9.2	3.8	25.6	26.2	45.5	2.6
sd	2.1	1.3	2.0	0.8	6.2	4.4	4.6	0.0

4 Agreement between models

4.1 Summary dataframe

4.1.1 Difference with AI ref for all SSP

cmip6 <- read.table("cmip6_AI.txt")

tab <- cmip6

list_ref <- list()

for(i in unique(tab$source)){
  tab.i <- filter(tab, source == i) 
  
  if(i == "CESM"){
    tab.Ai.ref <- filter(tab.i, period == "1970_2000") %>% pull(AI) %>% rep(8)
    tab.i$AI.ref <- tab.Ai.ref
    
    tab.pr.ref <- filter(tab.i, period == "1970_2000") %>% pull(spr) %>% rep(8)
    tab.i$pr.ref <- tab.pr.ref
    
    tab.tmean.ref <- filter(tab.i, period == "1970_2000") %>% pull(t) %>% rep(8)
    tab.i$t.ref <- tab.tmean.ref
    
  }else{
    tab.AI.ref <- filter(tab.i, period == "1970_2000") %>% pull(AI) %>% rep(9)
    tab.i$AI.ref <- tab.AI.ref
    
    tab.pr.ref <- filter(tab.i, period == "1970_2000") %>% pull(spr) %>% rep(9)
    tab.i$pr.ref <- tab.pr.ref
    
    tab.tmean.ref <- filter(tab.i, period == "1970_2000") %>% pull(t) %>% rep(9)
    tab.i$t.ref <- tab.tmean.ref
  }
  list_ref[[i]] <- tab.i  
}


cmip6ext <- bind_rows(list_ref, .id = "column_label") %>%
  mutate(diff.AI = AI - AI.ref, diff.pr = pr - pr.ref, diff.t = t - t.ref)
write.table(cmip6ext, "cmip6ext.txt")

cmip6 <- read.table("cmip6ext.txt")
cmip6s <- cmip6 %>% 
  group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, period, z) %>% 
  dplyr::summarise(t.mean = mean(t, na.rm = T), t.sd = sd(t, na.rm = T), # t in °C
                   pr.mean = mean(spr, na.rm = T) , pr.sd = sd(spr, na.rm = T),  # pr in mm/y
                   sfcWind.mean = mean(sfcWind, na.rm = T), sfcWind.sd = sd(sfcWind, na.rm = T),
                   Rn.mean = mean(Rn, na.rm = T), Rn.sd = sd(Rn, na.rm = T), # in MJ/m2/day
                   ea.mean = mean(ea, na.rm = T), ea.sd = sd(ea, na.rm = T), # in kPa
                   ET0 = mean(sET0, na.rm = T), ET0.sd = sd(sET0, na.rm = T), # in mm/y
                   AI.mean = mean(AI, na.rm = T), AI.sd = sd(AI, na.rm = T),
                   diff.AI.mean = mean(diff.AI, na.rm = T), diff.AI.sd = sd(diff.AI, na.rm = T),
                   hyperarid = sum(cat.AI == "Hyperarid"),
                   arid = sum(cat.AI == "Arid"),
                   semiarid = sum(cat.AI == "Semi-arid"),
                   drysubhumid = sum(cat.AI == "Dry subhumid") ,
                   humid = sum(cat.AI == "Humid"),
                   cold = sum(cat.AI == "Cold"),
                   na = sum(is.na(cat.AI)),
                   nb.models = n(),
                   nb.maj = max(hyperarid,arid,semiarid,drysubhumid,humid,cold,na),
                   prop.maj = nb.maj/nb.models,
                   plus.AI = table(diff.AI > 0)["TRUE"],
                   minus.AI = table(diff.AI < 0)["TRUE"]) %>% 
  ungroup() %>% as.data.frame()

cmip6s$cat.maj <- colnames(select(cmip6s, c("hyperarid","arid",'semiarid',"drysubhumid",'humid',"cold","na")))[max.col(select(cmip6s, c("hyperarid","arid",'semiarid',"drysubhumid",'humid',"cold","na")),ties.method = "first")]

write.table(cmip6s, "cmip6s.txt")

#test <- cmip6 %>% subset(lon == -10.5 & lat == 6.5 & Continent == "AFRICA" & model == "SSP245" & period == "2030_2060")

cmip6s <- read.table("cmip6s.txt")

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")

cmip6s$cat.AI <- cmip6s$AI.mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6s$cat.AI[which(cmip6s$ET0 < 400)] <- "Cold"

tab <- cmip6s
list_ref <- list()

for(i in unique(tab$lon)){
  tab.i <- filter(tab, lon == i)
  for(j in unique(tab$lat)){
    index <- paste(i,j,sep = ",")
    tab.ij <- filter(tab.i, lat == j)
    tab.ij$AIref <- filter(tab.ij, period == "1970_2000") %>% pull(AI.mean) %>% rep(9) #3 historical scenario + 2 x 3 future scenario
    tab.ij$catAIref <- filter(tab.ij, period == "1970_2000") %>% pull(cat.AI) %>% rep(9)
    tab.ij$catmajref <- filter(tab.ij, period == "1970_2000") %>% pull(cat.maj) %>% rep(9)
    list_ref[[index]] <- tab.ij 
  }
}

cmip6ss <- bind_rows(list_ref)

cmip6ss$anomalies <- with(cmip6ss, AI.mean - AIref)

cmip6ss$same.catAI <- with(cmip6ss, ifelse(cat.maj == cat.AI, "same", "diff"))

write.table(cmip6ss, "cmip6ss.txt")

cmip6s <- read.table("cmip6ss.txt")

# attributes a change from the studied period compared to the reference period 1970-2000. Read as: grid cell x is humid in 2030-2060 and Arid in 1970-2000, so it is WETTER

#cmip6s$change <-  with(cmip6s, paste(catAIref, cat.AI, collapse = " ", sep = " to "))

cmip6s$dryer <- with(cmip6s, 
          ifelse(catAIref == "Arid" & cat.AI %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
          ifelse(catAIref == "Arid" & cat.AI == "Hyperarid", "dryer",
          ifelse(catAIref == "Semi-arid" & cat.AI %in% c("Cold","Humid","Dry subhumid"), "wetter",
          ifelse(catAIref == "Semi-arid" & cat.AI %in% c("Hyperarid","Arid"), "dryer",
          ifelse(catAIref == "Dry subhumid" & cat.AI %in% c("Cold","Humid"), "wetter",
          ifelse(catAIref == "Dry subhumid" & cat.AI %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(catAIref == "Humid" & cat.AI == "Cold","wetter",
          ifelse(catAIref == "Humid" & cat.AI %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(catAIref == "Cold" & cat.AI %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(catAIref == "Hyperarid" & cat.AI %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
          ifelse(catAIref == cat.AI, "same",
          NA))))))))))))

cmip6s$catmajref <- revalue(cmip6s$catmajref, c("hyperarid" = "Hyperarid", "arid" = "Arid", "semiarid" = "Semi-arid", "drysubhumid" = "Dry subhumid", "humid" = "Humid", "cold" = "Cold"))

cmip6s$cat.maj <- revalue(cmip6s$cat.maj, c("hyperarid" = "Hyperarid", "arid" = "Arid", "semiarid" = "Semi-arid", "drysubhumid" = "Dry subhumid", "humid" = "Humid", "cold" = "Cold"))



cmip6s$dryer.maj <- with(cmip6s, 
        ifelse(catmajref == "Arid" & cat.maj %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
        ifelse(catmajref == "Arid" & cat.maj == "Hyperarid", "dryer",
        ifelse(catmajref == "Semi-arid" & cat.maj %in% c("Cold","Humid","Dry subhumid"), "wetter",
        ifelse(catmajref == "Semi-arid" & cat.maj %in% c("Hyperarid","Arid"), "dryer",
        ifelse(catmajref == "Dry subhumid" & cat.maj %in% c("Cold","Humid"), "wetter",
        ifelse(catmajref == "Dry subhumid" & cat.maj %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(catmajref == "Humid" & cat.maj == "Cold","wetter",
        ifelse(catmajref == "Humid" & cat.maj %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(catmajref == "Cold" & cat.maj %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(catmajref == "Hyperarid" & cat.maj %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
        ifelse(catmajref == cat.maj, "same",
        NA))))))))))))

4.2 Compare cat mean and cat maj

cat mean = category of aridity obtained by using the multimodel mean AI cat maj = category of aridity obtained by using the category majoritarily obtained in the 13 models

4.2.1 Where are cat mean and cat maj different

cmip6s %>% group_by(model, period) %>% summarise(maj.equal.mean = table(same.catAI)[2], maj.diff.mean = table(same.catAI)[1], prop.same = round(maj.equal.mean/n()*100,1)) %>% knitr::kable(col.names = c("Model","Period","Number of gridcells in which cat.maj = cat.mean", "Number of gridcells in which cat.maj /= cat.mean", "Proportion of cells with the same category (%)")) %>% kableExtra::kable_styling(bootstrap_options = "bordered")

Model	Period	Number of gridcells in which cat.maj = cat.mean	Number of gridcells in which cat.maj /= cat.mean	Proportion of cells with the same category (%)
SSP245	2030_2060	NA	21682	NA
SSP245	2070_2100	NA	21682	NA
SSP370	2030_2060	NA	21682	NA
SSP370	2070_2100	NA	21682	NA
SSP585	2030_2060	NA	21682	NA
SSP585	2070_2100	NA	21682	NA
historical	1850_1880	NA	21682	NA
historical	1970_2000	NA	21682	NA
historical	1985_2015	NA	21682	NA

cmip6s$cat.AI <- factor(cmip6s$cat.AI, levels =  c("Hyperarid","Arid","Semi-arid","Dry subhumid", "Humid","Cold"))
cmip6s$cat.maj <- factor(cmip6s$cat.maj, levels =  c("Hyperarid","Arid","Semi-arid","Dry subhumid", "Humid","Cold"))

ggplot(subset(cmip6s, cat.AI != "Cold"))+
  geom_col(aes(x= period, y = cat.AI, fill = cat.AI), just = -0.2, width = 0.3)+
  geom_col(aes(x = period, y = cat.maj, fill = cat.maj), just = 1.2, width = 0.3)+
  scale_fill_manual(values = col.cat)+
  labs(fill = "", y ="")+
  facet_grid(rows = vars(model), switch = "y")+
  theme_minimal()+theme(axis.text.y = element_blank())

tab <- rbind(table(cmip6s$cat.AI)/length(cmip6s$cat.AI)*100, table(cmip6s$cat.maj)/length(cmip6s$cat.maj)*100) %>% t() %>% as.data.frame() %>% setNames(c("Cat mean", "Cat maj"))
knitr::kable(tab, digits =2, caption = "Proportion of gridcells in each aridity category") %>% kable_styling(bootstrap_options = "bordered")

Proportion of gridcells in each aridity category

	Cat mean	Cat maj
Hyperarid	3.51	4.18
Arid	6.25	7.43
Semi-arid	7.63	9.19
Dry subhumid	3.39	0.84
Humid	28.75	28.53
Cold	47.90	47.23

map_list <- list()

for(i in unique(cmip6s$model)){
  df <- subset(cmip6s, model == i)
  for(j in unique(df$period)){
    index <- paste(i, j, sep = "")
    g <- ggplot()+geom_raster(data = subset(df, period == j & same.catAI == "diff"), aes(x = lon, y = lat, fill = cat.AI))+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste(i, j , sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }
}

col1 <- ggarrange(plotlist = map_list, nrow = 9, ncol = 1, common.legend = T, legend = "bottom") %>% annotate_figure(top = "cat.mean for gridcells in which\ncat.maj and cat.mean are different")

map_list2 <- list()

for(i in unique(cmip6s$model)){
  df <- subset(cmip6s, model == i)
  for(j in unique(df$period)){
    index <- paste(i, j, sep = "")
    g <- ggplot()+geom_raster(data = subset(df, period == j & same.catAI == "diff"), aes(x = lon, y = lat, fill = cat.maj))+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste(i, j , sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list2[[index]] <- g
  }
}

col2 <- ggarrange(plotlist = map_list2, nrow = 9, ncol = 1, common.legend = T, legend = "bottom") %>% annotate_figure(top = "cat.maj for gridcells in which\ncat.maj and cat.mean are different")

ggarrange(col1,col2,ncol = 2, common.legend = T)

4.2.2 Plot cat maj and cat mean for reference period

ggpubr::ggarrange(plotlist = list(
  
  ggplot() + 
    geom_raster(data = subset(cmip6s, period == "1970_2000"), aes(x=lon, y = lat,  fill = cat.maj))+
    geom_point(data = subset(cmip6s, period == "1970_2000" & nb.maj > 10), aes(x =lon, y = lat), shape = "'", size = 1, col = "grey20")+
    scale_fill_manual(values = col.cat)+
    labs(title = "Category chosen in majority by models", fill = "")+
    ylim(-55,90)+
    theme_void(base_size = 15)+
    theme(legend.position = "none"),
  
  ggplot() + 
    geom_raster(data = subset(cmip6s, period == "1970_2000"), aes(x=lon, y = lat,  fill = cat.AI))+
    geom_point(data = subset(cmip6s, period == "1970_2000" & nb.maj > 10), aes(x =lon, y = lat), shape = "'", size = 1, col = "grey20")+
    scale_fill_manual(values = col.cat)+
    labs(title = "Multimodel average categories of AI", fill = "")+
    ylim(-55,90)+
    theme_void(base_size = 15)+
    theme(legend.position = "none")),
  
  ncol = 2)

map_list <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.maj))+
    geom_point(data = subset(cmip6s, period == i & model == "historical" & nb.maj > 10), aes(x= lon, y = lat), shape = "'", col = "black")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
}


map_list2 <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI))+
    geom_point(data = subset(cmip6s, period == i & model == "historical" & nb.maj > 10), aes(x= lon, y = lat), shape = "'", col = "black")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list2[[i]] <- g
}


ggarrange(
  annotate_figure(ggarrange(plotlist = map_list, nrow = 3, common.legend = T, legend = "bottom"), top = "Majority category among models"),
  annotate_figure(ggarrange(plotlist = map_list2, nrow = 3, common.legend = T, legend = "bottom"), top = "Multimodel average category"),
  ncol = 2, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.maj))+
      #   geom_point(data = subset(cmip6s, period == i & model == "historical" & change.cat == 1 & nb.maj > 10), aes(x= lon, y = lat), shape = "'")+
      scale_fill_manual(values = col.cat, na.translate = F)+
      labs(fill = "Cat AI maj", title = paste(i, j, sep = ", "))+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}


map_list2 <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.AI))+
      #   geom_point(data = subset(cmip6s, period == i & model == "historical" & change.cat == 1 & nb.maj > 10), aes(x= lon, y = lat), shape = "'")+
      scale_fill_manual(values = col.cat, na.translate = F)+
      labs(fill = "Cat AI mean", title = paste(i, j, sep = ", "))+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list2[[index]] <- g
  }}

ggarrange(plotlist = map_list, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

ggarrange(plotlist = map_list2, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

4.3 Multimodel maps of cat mean

cmip6s$cat.AI <- factor(cmip6s$cat.AI, levels = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))

map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + 
    geom_tile(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI))+
    geom_point(data = subset(cmip6s,nb.maj > 10), aes(x = lon, y = lat), shape = "'", col = "grey60")+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = paste("Category AI",i , sep = " "), fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + 
      geom_tile(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.AI))+
      geom_point(data = subset(cmip6s,nb.maj > 10), aes(x = lon, y = lat), shape = ".", col = "grey20", size = 0.1)+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste("Category AI",index, sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

5 Maps multimodel AI

quantile(cmip6s$AI.mean, probs = seq(0,1,0.1), na.rm = T) 
##           0%          10%          20%          30%          40%          50% 
## 2.902197e-03 1.993698e-01 5.922059e-01 1.019621e+00 1.457766e+00 2.127218e+00 
##          60%          70%          80%          90%         100% 
## 3.578980e+00 1.168470e+01 2.438867e+01 4.856852e+01 5.329875e+06

ai.breaks <- c(0,0.03,0.2,0.5,0.65,1,10,20,30)
colscale <- colorvec[-c(5,6)]

map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = AI.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = ai.breaks, palette = function(x) colscale,
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Aridity index")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = AI.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = ai.breaks, palette = function(x) colscale,
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Aridity index", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

6 AI changes over time and scenarios

v<-NULL
v2<-NULL
for (N in 5:30) {
  v<-c(v,sum(sapply(floor(2*N/3):N,function(i) choose(N,i)*(1/2)^N*2)))# proba d'avoir 2/3 de même signe à partir de N modèles, lorsque N change
  v2<-c(v2,sum(sapply(floor(2*N/3):N,function(i) choose(N,i)*(1/2)^N)))# proba d'avoir 2/3 de signe >0 à partir de N modèles, lorsque N change
}


plot(5:30,v,type="l",ylab="Proba(>2/3 models...)",xlab="Nb of models")
lines(5:30,v2,col="red")

7 Changes AI for summary CMIP6

7.1 Boxplot of changes AI

ggpubr::ggarrange(plotlist = list(
  ggplot(subset(cmip6s, anomalies < 1 & anomalies > -1 & model == "historical"))+
    geom_boxplot(aes(x=period, y = - anomalies, col = period), outlier.size = 1)+
    scale_color_manual(values = c("#1A9850","#66BD63","#A6D96A"))+
    labs(x = "", y = "difference of AI between each period and\nthe reference period 1970-2000")+
    facet_grid(cols = vars(model))+
    theme_minimal() + theme(legend.position = "none"),
  ggplot(subset(cmip6s, anomalies < 1 & anomalies > -1 & model != "historical"))+
    geom_boxplot(aes(x=period, y = anomalies, col = period), outlier.size = 1)+
    scale_color_manual(values = c("#FDAE61","#F46D43"))+
    scale_y_continuous(position = "right")+
    labs(y = "", x = "")+
    facet_grid(cols = vars(model))+
    theme_minimal() + theme(legend.position = "none")),
  widths = c(1.5,3),
  ncol = 2, nrow = 1
)

7.2 Strictly wetter or dryer (reference period: 1970-2000)

test <- cmip6s %>% select(c("lon","lat","model","period","dryer","cat.AI","catAIref","AI.mean", "AIref", "anomalies","plus.AI", "minus.AI"))

View(subset(test, dryer == "dryer"))

ggplot(subset(cmip6s, lat > -55 & !is.na(dryer))) + 
  geom_boxplot(aes(y = cat.AI, x = diff.AI.mean, col = catAIref), outliers = F) +
  scale_color_manual(values = col.cat)+
  facet_grid(cols = vars(dryer))+
  labs(y = "Aridity category", x = "Difference between AI and AI ref", col = "Reference aridity category")+
  theme_bw()

ggplot(subset(cmip6s, lat > -55 & !is.na(dryer))) + 
  geom_boxplot(aes(y = cat.AI, x = anomalies, col = catAIref), outliers = F) +
  scale_color_manual(values = col.cat)+
  facet_grid(cols = vars(dryer))+
  labs(y = "Aridity category", x = "Difference between AI and AI ref", col = "Reference aridity category")+
  theme_bw()

ggplot(cmip6s)+
  geom_point(data = subset(cmip6s, plus.AI > 10), aes(x = lon, y = lat, col = "plus"), shape = "+")+
  geom_point(data = subset(cmip6s, minus.AI > 10), aes(x = lon, y = lat, col = "minus"), shape = "-")+
  theme_minimal()

filter(cmip6s, plus.AI > 10 & minus.AI > 10) %>% View()

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
  index <- paste(i, j, sep = ", ")
  g <- ggplot()+
  geom_point(data = subset(cmip6s, model == i & period == j & plus.AI > 10), aes(x = lon, y = lat, col = "plus"), shape = "+")+
  geom_point(data = subset(cmip6s, model == i & period == j & minus.AI > 10), aes(x = lon, y = lat, col = "minus"), shape = "-")+
    scale_color_manual(values = colorvec[c(1,11)])+
    labs(color = "", x = "", y = "", title = index)+
  theme_minimal()
  map_list[[index]] <- g
  }
}

ggarrange(plotlist = map_list, ncol = 2, nrow = 3, common.legend = T)

map_list <- list()

for(i in c("1850_1880","1985_2015")){
    index <- paste(i, j, sep = " , ")
  g <- ggplot() + geom_tile(data = subset(cmip6s, period == i & model == "historical" & anomalies > 0), aes(x=lon, y = lat,  fill = "wetter"))+
    geom_tile(data = subset(cmip6s, period == i & model == "historical" & anomalies < 0), aes(x=lon, y = lat,  fill = "dryer"))+
    geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = ".")+
    borders(colour = "grey60")+
    scale_fill_manual(values =  c("#F46D43", "#74ADD1"), na.translate = F)+
    labs(title = paste("AI changes between",i , "and 1970-2000", sep = " "), fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
  }

ggarrange(plotlist = map_list, ncol = 2)

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & anomalies > 0), aes(x=lon, y = lat,  fill = "wetter"))+
      geom_tile(data = subset(cmip6s, period == j & model == i & anomalies < 0), aes(x=lon, y = lat,  fill = "dryer"))+
       geom_point(data = subset(cmip6s, period == j & model == i & anomalies > 0 & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey40", size = 0.1)+
       geom_point(data = subset(cmip6s, period == j & model == i & anomalies < 0 & minus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey40", size = 0.1)+
      borders(colour = "grey60")+
      scale_fill_manual(values = c("#F46D43", "#74ADD1"), na.translate = F)+
      labs(fill = "Compared to 1970-2000", title = i)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

7.3 Wetter or dryer, binned (reference period: 1970-2000)

ano.quantiles <- cmip6s %>% group_by(model, period) %>% summarise(q10 = quantile(anomalies, probs = 0.1, na.rm =T), 
                                                                  q25 = quantile(anomalies, probs = 0.25, na.rm =T), 
                                                                  q50 = quantile(anomalies, probs = 0.5, na.rm =T), 
                                                                  q75 = quantile(anomalies, probs = 0.75, na.rm =T), 
                                                                  q90 = quantile(anomalies, probs = 0.9, na.rm =T)) 

ano.quantiles <- cmip6s %>% summarise(q10 = quantile(anomalies, probs = 0.1, na.rm =T), 
                                      q25 = quantile(anomalies, probs = 0.25, na.rm =T), 
                                      q50 = quantile(anomalies, probs = 0.5, na.rm =T), 
                                      q75 = quantile(anomalies, probs = 0.75, na.rm =T), 
                                      q90 = quantile(anomalies, probs = 0.9, na.rm =T)) 
qbreaks <- c(-5, -1,-0.1, -0.01, 0, 0.01, 0.1, 1,5)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")

“Anomalies” is the multimodel AI average minus the reference AI (period 1970-2000). More blue: AI mean is superior to AI ref, hence wetter than the ref. More red: the anomalie is negative, hence AI mean is inferior to AI ref, hence dryer.

AS OF NOW THE SAME QUANTILES ARE USED FOR ALL MAPS.

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  #quant = subset(ano.quantiles, model == "historical" & period == i)
  #b = quant[1,c(3:7)] %>% as.numeric
  b = qbreaks
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = anomalies))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    b = qbreaks
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = anomalies))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Compared to 1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

7.4 Wetter or dryer, binned, % (reference period: 1970-2000)

cmip6s$anomalies.percent <- with(cmip6s, (AI.mean-AIref)/AIref*100)


pbreaks <- c(-40, -30,-20, -10, 0, 10, 20, 30, 40)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")

“Anomalies” is the multimodel AI average minus the reference AI (period 1970-2000). More blue: AI mean is superior to AI ref, hence wetter than the ref. More red: the anomalie is negative, hence AI mean is inferior to AI ref, hence dryer.

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  #quant = subset(ano.quantiles, model == "historical" & period == i)
  #b = quant[1,c(3:7)] %>% as.numeric
  b = pbreaks
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = anomalies.percent))+
    borders(colour = "grey60")+
           geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey40", size = 0.1)+
       geom_point(data = subset(cmip6s,  period == i & model == "historical"  & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey40", size = 0.1)+
    binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 2, common.legend = T, legend = "bottom") %>% annotate_figure(top = "% change of AI index compared to 1970-2000")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    b = pbreaks
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = anomalies.percent))+
      borders(colour = "grey60", size = 0.5)+
                 geom_point(data = subset(cmip6s, period == j & model == i & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey40", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == j & model == i & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey40", size = 0.5)+
      binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "% change of AI index\ncompared to 1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}


ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

8 Difference between future and reference category of AI, using cat maj

8.1 Detailed changes and colors cat maj

cmip6s$change.catmaj <- with(cmip6s, ifelse(cat.maj == catmajref, 0, 1))

cmip6s$change.cat.maj <- with(cmip6s, paste(catmajref, cat.maj, sep = " to "))

cmip6s$change.cat.maj <- with(cmip6s, ifelse(catmajref == "Arid" & cat.maj %in% c("Cold","Humid"),"Arid to cold or humid",
         ifelse(catmajref == "Arid" & cat.maj %in% c("Semi-Arid","Dry subhumid"), "Arid to semi-arid or dry subhumid",
         ifelse(catmajref == "Arid" & cat.maj == "Hyperarid", "Arid to hyperarid",
         ifelse(catmajref == "Semi-arid" & cat.maj == "Dry subhumid", "Semi-arid to dry subhumid",
         ifelse(catmajref == "Semi-arid" & cat.maj == "Arid", "Semi-arid to arid",
         ifelse(catmajref == "Semi-arid" & cat.maj == "Hyperarid", "Semi-arid to hyperarid",
         ifelse(catmajref == "Semi-arid" & cat.maj %in% c("Humid","Cold"), "Semi-arid to humid or cold",
         ifelse(catmajref == "Dry subhumid" & cat.maj %in% c("Semi-arid","Arid","Hyperarid"), "Dry subhumid to arid, semi-arid or hyperarid",
         ifelse(catmajref == "Dry subhumid" & cat.maj %in% c("Humid","Cold"), "Dry subhumid to humid or cold",
         ifelse(catmajref == "Humid" & cat.maj == "Cold","Humid to cold",
         ifelse(catmajref == "Humid" & cat.maj %in% c("Dry subhumid", "Semi-arid"), "Humid to dry-subhumid or semi-arid",
         ifelse(catmajref == "Humid" & cat.maj %in% c("Arid", "Hyperarid"), "Humid to arid or hyperarid",
         ifelse(catmajref == "Cold" & cat.maj == "Humid", "Cold to humid", 
         ifelse(catmajref == "Cold" & cat.maj %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "Cold to dryland",
         ifelse(catmajref == "Hyperarid" & cat.maj == "Arid", "Hyperarid to arid",
         ifelse(catmajref == "Hyperarid" & cat.maj %in% c("Semi-arid","Dry-subhumid"),"Hyperarid to semi arid or dry subhumid",
         ifelse(catmajref == "Hyperarid" & cat.maj %in% c("Humid", "Cold"), "Hyperarid to humid or cold", NA))))))))))))))))))

cmip6s$change.cat.maj <- factor(cmip6s$change.cat.maj, 
                                levels = c("Arid to hyperarid","Semi-arid to arid", "Dry subhumid to arid, semi-arid or hyperarid", "Humid to dry-subhumid or semi-arid", "Cold to humid","Hyperarid to arid", "Arid to semi-arid or dry subhumid", "Dry subhumid to humid or cold", "Semi-arid to dry subhumid", "Semi-arid to humid or cold"))
table(cmip6s$change.cat.maj)
## 
##                            Arid to hyperarid 
##                                          491 
##                            Semi-arid to arid 
##                                         1462 
## Dry subhumid to arid, semi-arid or hyperarid 
##                                          649 
##           Humid to dry-subhumid or semi-arid 
##                                         3088 
##                                Cold to humid 
##                                         3700 
##                            Hyperarid to arid 
##                                          301 
##            Arid to semi-arid or dry subhumid 
##                                            4 
##                Dry subhumid to humid or cold 
##                                          201 
##                    Semi-arid to dry subhumid 
##                                          236 
##                   Semi-arid to humid or cold 
##                                          381

colors_changes <- c("Hyperarid to arid" = "#E7D7B6", 
                    "Arid to semi-arid or dry subhumid" = colorvec[7],
                    "Dry subhumid to humid or cold" = colorvec[8],
                    "Semi-arid to dry subhumid" = colorvec[9],
                    "Semi-arid to humid or cold" = colorvec[10],
                    "Humid to cold" = colorvec[11],
                    "Cold to humid" = "#E7D7B6", 
                    "Humid to dry-subhumid or semi-arid" = colorvec[5],
                    "Dry subhumid to arid, semi-arid or hyperarid" = colorvec[4],
                    "Semi-arid to arid" = colorvec[3],
                    "Arid to hyperarid" = colorvec[1])

8.2 Changes towards dryer categories cat maj

8.2.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + 
    geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & dryer.maj == "dryer"), 
                aes(x=lon, y = lat,  fill = change.cat.maj))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}


ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

8.2.1.1 Future

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmaj == 1 & dryer.maj == "dryer"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
      #   geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & nb.maj > 10), aes(x= lon, y = lat, fill = "'"))+
      borders(colour = "grey60")+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      #  scale_fill_viridis_d()+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

8.3 Changes towards wetter categories cat maj

8.3.0.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & dryer.maj == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
    borders(colour = "grey60", ylim = c(-40,40))+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+
    ylim(-50,50)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

8.3.0.2 Future

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmaj == 1 & dryer.maj == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
      borders(colour = "grey60", ylim = c(-40,40))+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-50,50)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

9 Difference between future and reference category of AI, using cat mean

9.1 Detailed changes and colors

cmip6s$change.catmean <- with(cmip6s, ifelse(cat.AI == catAIref, 0, 1))

cmip6s$change.cat.mean <- with(cmip6s, paste(catAIref, cat.AI, sep = " to "))
unique(cmip6s$change.cat)
## NULL

# cmip6s$change.cat.mean <- with(cmip6s, 
#       ifelse(catAIref == "Arid" & cat.AI %in% c("Cold","Humid"),"Arid to cold or humid",
#       ifelse(catAIref == "Arid" & cat.AI %in% c("Semi-Arid","Dry subhumid"), "Arid to semi-arid or dry subhumid",
#       ifelse(catAIref == "Arid" & cat.AI == "Hyperarid", "Arid to hyperarid",
#       ifelse(catAIref == "Semi-arid" & cat.AI == "Dry subhumid", "Semi-arid to dry subhumid",
#       ifelse(catAIref == "Semi-arid" & cat.AI == "Arid", "Semi-arid to arid",
#       ifelse(catAIref == "Semi-arid" & cat.AI == "Hyperarid", "Semi-arid to hyperarid",
#       ifelse(catAIref == "Semi-arid" & cat.AI %in% c("Humid","Cold"), "Semi-arid to humid or cold",
#       ifelse(catAIref == "Dry subhumid" & cat.AI %in% c("Semi-arid","Arid","Hyperarid"), "Dry subhumid to arid, semi-arid or hyperarid",
#       ifelse(catAIref == "Dry subhumid" & cat.AI %in% c("Humid","Cold"), "Dry subhumid to humid or cold",
#       ifelse(catAIref == "Humid" & cat.AI == "Cold","Humid to cold",
#       ifelse(catAIref == "Humid" & cat.AI %in% c("Dry subhumid","Semi-arid"), "Humid to dry-subhumid or semi-arid",
#       ifelse(catAIref == "Humid" & cat.AI %in% c("Arid", "Hyperarid"), "Humid to arid or hyperarid",
#       ifelse(catAIref == "Cold" & cat.AI == "Humid", "Cold to humid", 
#       ifelse(catAIref == "Cold" & cat.AI %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "Cold to dryland",
#       ifelse(catAIref == "Hyperarid" & cat.AI == "Arid", "Hyperarid to arid",
#       ifelse(catAIref == "Hyperarid" & cat.AI %in% c("Semi-arid","Dry-subhumid"),"Hyperarid to semi arid or dry subhumid",
#       ifelse(catAIref == "Hyperarid" & cat.AI %in% c("Humid", "Cold"), "Hyperarid to humid or cold", 
#       NA))))))))))))))))))

table(cmip6s$change.cat.mean)
## 
##                 Arid to Arid            Arid to Hyperarid 
##                        11085                          638 
##            Arid to Semi-arid                 Cold to Cold 
##                          346                        95818 
##                Cold to Humid Dry subhumid to Dry subhumid 
##                         4370                         4059 
##        Dry subhumid to Humid    Dry subhumid to Semi-arid 
##                          297                         2214 
##                Humid to Cold        Humid to Dry subhumid 
##                          148                         2448 
##               Humid to Humid           Humid to Semi-arid 
##                        52938                          176 
##            Hyperarid to Arid       Hyperarid to Hyperarid 
##                          248                         6403 
##                     NA to NA            Semi-arid to Arid 
##                         5130                         1196 
##    Semi-arid to Dry subhumid           Semi-arid to Humid 
##                          276                            1 
##       Semi-arid to Semi-arid 
##                        12558
cmip6s$change.cat.mean <- factor(cmip6s$change.cat.mean, 
                                  levels = c("Arid to Hyperarid","Semi-arid to Arid", "Dry subhumid to Semi-arid", "Humid to Semi-arid", "Humid to Dry subhumid", "Cold to Humid","Hyperarid to arid", "Arid to Semi-arid", "Dry subhumid to Humid", "Semi-arid to Dry subhumid", "Semi-arid to Humid"))

colors_changes <- c("Hyperarid to Arid" = "#E7D7B6", 
                     "Arid to Semi-arid" = colorvec[7],
                     "Dry subhumid to Humid" = colorvec[8],
                     "Semi-arid to Dry subhumid" = colorvec[9],
                     "Semi-arid to Humid" = colorvec[10],
                     "Humid to Cold" = colorvec[11],
                     "Cold to Humid" = "#E7D7B6",
                     "Humid to Dry subhumid" = colorvec[5],
                     "Humid to Semi-arid" = colorvec[4],
                     "Dry subhumid to Semi-arid" = colorvec[3],
                     "Semi-arid to Arid" = colorvec[2],
                     "Arid to Hyperarid" = colorvec[1])

9.2 Changes towards dryer categories cat mean

9.2.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + 
    geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmean == 1 & dryer == "dryer"), 
                aes(x=lon, y = lat,  fill = change.cat.mean))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

9.2.1.1 Future

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmean == 1 & dryer == "dryer"), 
                              aes(x=lon, y = lat,  fill = change.cat.mean))+
      borders(colour = "grey60")+
      scale_fill_manual(values = colors_changes)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggarrange(plotlist = map_list, ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

9.3 Changes towards wetter categories cat mean

9.3.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + geom_tile(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & dryer == "wetter"), 
                            aes(x=lon, y = lat,  fill = change.cat.mean))+
    borders(colour = "grey60", ylim = c(-40,40))+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-50,50)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

9.3.2 Future

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmean == 1 & dryer == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.mean))+
      borders(colour = "grey60", ylim = c(-40,40))+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

#s Tables

9.4 Changes of aridity category using cat mean

tab <- cmip6s %>% group_by(model, period, Continent) %>% summarise(Cold = table(cat.AI)[1], Hyperarid = table(cat.AI)[2], Arid = table(cat.AI)[3], Semiarid = table(cat.AI)[4], Drysubhumid = table(cat.AI)[5], Humid = table(cat.AI)[6], total = n()) 
write.table(tab, "tab.catAI.txt")

tab.future <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(period + model ~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

kable(tab.future[order(tab.future$model, decreasing = F),]) %>% kable_styling(bootstrap_options = "bordered") %>%
  column_spec(c(7,10), italic = T, include_thead = T) %>% row_spec(2, bold = T) 

	model	period	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
1	historical	1850_1880	3.4	6.0	7.1	3.3	19.8	28.8	48.9	2.6
2	historical	1970_2000	3.3	6.0	7.0	3.3	19.6	27.8	49.9	2.6
3	historical	1985_2015	3.3	6.1	6.9	3.3	19.6	27.3	50.6	2.6
4	SSP245	2030_2060	3.5	6.3	7.8	3.3	20.9	28.8	47.7	2.6
7	SSP245	2070_2100	3.6	6.4	7.8	3.5	21.3	29.4	46.6	2.6
5	SSP370	2030_2060	3.2	6.7	7.6	3.0	20.5	28.8	48.1	2.6
8	SSP370	2070_2100	4.0	6.3	8.3	3.7	22.3	28.5	46.6	2.6
6	SSP585	2030_2060	3.6	6.2	7.9	3.5	21.2	29.3	47.0	2.6
9	SSP585	2070_2100	3.8	6.4	8.4	3.6	22.2	30.1	45.1	2.6

# ggplot(tab.future, aes(x = period))+
#   geom_point(aes(y = Hyperarid, col = "Hyperarid", shape = model))+
#   geom_line(aes(y=Hyperarid, col = "Hyperarid", group = model, lty = model))+
#   geom_point(aes(y = Arid, col = "Arid", shape = model))+
#   geom_line(aes(y=Arid, col = "Arid", group = model, lty = model))+
#   geom_point(aes(y = get('Semi-arid'), col = "Semi-arid", shape = model))+
#   geom_line(aes(y=get('Semi-arid'), col = "Semi-arid", group = model, lty = model))+
#   geom_point(aes(y = get('Dry subhumid'), col = "Dry subhumid", shape = model))+
#   geom_line(aes(y=get('Dry subhumid'), col = "Dry subhumid", group = model, lty = model))+
#   geom_point(aes(y = Humid, col = "Humid", shape = model))+
#   geom_line(aes(y= Humid, col = "Humid", group = model, lty = model))+
#   scale_color_manual(values = col.cat)+
#   theme_minimal()

cmip6$cat.AI <- factor(cmip6$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))
tab.percent <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","POLAR","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() %>% 
  group_by(period, model, cat.AI) %>% 
  summarise(mmmean = mean(count, na.rm = T), mmsd = sd(count, na.rm = T)) %>%
  mutate(percent = round(mmmean/sum(mmmean)*100, 1), sd.percent = round(mmsd/sum(mmmean)*100,1)) 
write.table(tab.percent, "tab.percent.txt")

tab.percent <- read.table("tab.percent.txt")

df245 <- rbind(subset(tab.percent, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.percent, model == "SSP245")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df245$cat.AI <- factor(df245$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 2-4.5", y = "%", x = "", col = "", fill = "")+
  theme_minimal()

df370 <- rbind(subset(tab.percent, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.percent, model == "SSP370")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df370$cat.AI <- factor(df370$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 3-7.0", y = "%", x = "", col = "", fill = "")+
  theme_minimal()


df585 <- rbind(subset(tab.percent, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.percent, model == "SSP585")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df585$cat.AI <- factor(df585$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 5-8.5", y = "%", x = "", col = "", fill = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

9.5 Changes of aridity category using cat maj

tab.maj <- cmip6s %>% group_by(model, period, Continent) %>% summarise(Cold = table(cat.maj)[1], Hyperarid = table(cat.maj)[2], Arid = table(cat.maj)[3], Semiarid = table(cat.maj)[4], Drysubhumid = table(cat.maj)[5], Humid = table(cat.maj)[6], total = n()) 
write.table(tab.maj, "tab.catmaj.txt")

tab.future.maj <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, cat.maj) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(period + model ~cat.maj) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

kable(tab.future.maj[order(tab.future.maj$model, decreasing = F),]) %>% kable_styling(bootstrap_options = "bordered") %>%
  column_spec(c(7,10), italic = T, include_thead = T) %>% row_spec(2, bold = T) 

	model	period	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
1	historical	1850_1880	4.2	7.1	8.7	0.8	20.8	28.3	48.3	2.6
2	historical	1970_2000	4.1	7.1	8.6	0.7	20.5	28.1	48.9	2.6
3	historical	1985_2015	4.0	7.0	8.4	0.9	20.3	27.3	49.9	2.6
4	SSP245	2030_2060	4.0	7.8	9.3	0.8	21.9	28.4	47.2	2.6
7	SSP245	2070_2100	4.2	7.7	9.4	0.9	22.2	29.1	46.1	2.6
5	SSP370	2030_2060	4.1	7.2	9.4	0.8	21.5	28.6	47.3	2.6
8	SSP370	2070_2100	4.5	7.7	10.0	0.9	23.1	28.4	45.9	2.6
6	SSP585	2030_2060	4.2	7.5	9.6	1.0	22.3	28.8	46.3	2.6
9	SSP585	2070_2100	4.5	7.9	9.4	0.9	22.7	30.0	44.6	2.6

tab.flow.catmaj <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","POLAR","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, cat.maj) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL)

df245 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow.catmaj, model == "SSP245")) %>%
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 2-4.5", y = "%", x = "")+
  theme_minimal()

df370 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow.catmaj, model == "SSP370")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 3-7.0", y = "%", x = "")+
  theme_minimal()


df585 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow.catmaj, model == "SSP585")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 5-8.5", y = "%", x = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

9.6 Changes by continent, cat mean

cmip6$cat.AI <- factor(cmip6$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

tab.percent.cont <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","POLAR","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, Continent, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() %>% 
  group_by(period, model, Continent, cat.AI) %>% 
  summarise(mmmean = mean(count, na.rm = T), mmsd = sd(count, na.rm = T)) %>%
  mutate(percent = round(mmmean/sum(mmmean)*100, 1), sd.percent = round(mmsd/sum(mmmean)*100,1)) 
write.table(tab.percent.cont, "tab.percent.continent.txt")

tab.cont <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","POLAR","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, Continent, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() 

list_ref <- list()

for(i in unique(tab.cont$source)){
  tab.cont.i <- filter(tab.cont, source == i)
  tab.cont.i.ref <- filter(tab.cont.i, period == "1970_2000")
  
  for(j in unique(tab.cont.i$Continent)){
    tab.cont.i.j <- filter(tab.cont.i, Continent == j)
    tab.cont.i.j.ref <- filter(tab.cont.i.ref, Continent == j)
    
    for(k in unique(tab.cont.i$cat.AI)){
      
      index <- paste(i, j, k, sep = "_")
      tab.cont.i.j.k <- filter(tab.cont.i.j, cat.AI == k)
      
      lines <- ifelse(dim(tab.cont.i.j.k)[1] > 0, dim(tab.cont.i.j.k)[1] - 1, 1)
      
      tab.cont.i.j.k.ref <- subset(tab.cont.i.j.ref, cat.AI == k)
      
      if(dim(tab.cont.i.j.k.ref)[1] == 0){
        names1 <- names(tab.cont.i.j.k.ref)
        line1 <- c("1970_2000","historical", i, j, k, 0)
        tab.cont.i.j.k.ref <- rbind(tab.cont.i.j.k.ref, line1) %>%
          setNames(names1)
        
      }else{
        tab.cont.i.j.k.ref <- tab.cont.i.j.k.ref
      }
      
      tab.cont.i.j.k.ref <- tab.cont.i.j.k.ref %>% rbind(tab.cont.i.j.k.ref[rep(1,lines),])
      tab.cont.i.j.k$diff <- tab.cont.i.j.k$count - as.numeric(tab.cont.i.j.k.ref$count)  
      
      list_ref[[index]] <- tab.cont.i.j.k
    }
  }
}


tab.cont.diff <- bind_rows(list_ref, .id = "column_label") %>%
  group_by(period, model, Continent, cat.AI) %>% 
  summarise(diff.mean = mean(diff, na.rm = T), diff.sd = sd(diff, na.rm = T))
write.table(tab.cont.diff, "tab.cont.diff.txt")

tab.percent.cont <- read.table("tab.percent.continent.txt")

tab_list_percent <- list()
tab_list_sd <- list()

for(i in unique(tab.percent.cont$Continent)){
  tab.i <- subset(tab.percent.cont, Continent == i)
  
  df.percent <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "percent")
  
  df.sd <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "sd.percent")
  # sd for a sum: square root of sum of squared sd
  
  drycats <- which(names(df.percent) %in% c("Hyperarid","Arid","Semi-arid","Dry subhumid"))
  
  df.percent <- df.percent %>% mutate("Sum drylands" = rowSums(.[drycats], na.rm = T))
  df.sd <- df.sd %>% mutate("Sum drylands" = round(sqrt(rowSums((.[drycats])^2, na.rm = T)),2))
  
  missing <- setdiff(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"), names(df.percent))
  
  df.percent[, missing] <- 0
  df.sd[,missing] <- 0
  
  df.percent <- df.percent %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.percent[is.na(df.percent)] <- 0
  
  tab_list_percent[[i]] <- df.percent
  
  df.sd <- df.sd %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.sd[is.na(df.sd)] <- 0
  
  tab_list_sd[[i]] <- df.sd
}

k_list <- list()
kdiff_list <- list()

for(i in names(tab_list_percent)){
  
  df <- tab_list_percent[[i]] %>% select(c("model", "period"))
  
  prop.ref <- tab_list_percent[[i]] %>% filter(period == "1970_2000")
  prop.ref <- prop.ref %>% rbind(prop.ref[rep(1,8),])
  df2 <- tab_list_percent[[i]][,c(3:10)] - prop.ref[,c(3:10)] %>% round(1)
  diff.prop.percent <- cbind(df, df2)
  
  sd.ref <- tab_list_sd[[i]] %>% filter(period == "1970_2000")
  sd.ref <- sd.ref %>% rbind(sd.ref[rep(1,8),])
  df3 <- sqrt(tab_list_sd[[i]][,c(3:10)]^2 + sd.ref[,c(3:10)]^2) %>% round(1)
  diff.sd.percent <- cbind(df, df3)
  
  df$Hyperarid <- paste(tab_list_percent[[i]]$Hyperarid, tab_list_sd[[i]]$Hyperarid, sep = " ± ")
  df$Arid <- paste(tab_list_percent[[i]]$Arid, tab_list_sd[[i]]$Arid, sep = " ± ")
  df$'Semi-Arid' <- paste(tab_list_percent[[i]][,5], tab_list_sd[[i]][,5], sep = " ± ")
  df$'Dry subhumid' <- paste(tab_list_percent[[i]][,6], tab_list_sd[[i]][,6], sep = " ± ")
  df$'Sum drylands' <- paste(tab_list_percent[[i]][,7], tab_list_sd[[i]][,7], sep = " ± ")
  df$Humid <- paste(tab_list_percent[[i]]$Humid, tab_list_sd[[i]]$Humid, sep = " ± ")
  df$Cold <- paste(tab_list_percent[[i]]$Cold, tab_list_sd[[i]]$Cold, sep = " ± ")
  
  df.diff <- df
  df.diff$Hyperarid <- paste(round(diff.prop.percent$Hyperarid, 1), round(diff.sd.percent$Hyperarid, 1), sep = " ± ")
  df.diff$Arid <- paste(round(diff.prop.percent$Arid, 1), round(diff.sd.percent$Arid, 1), sep = " ± ")
  df.diff$'Semi-Arid' <- paste(round(diff.prop.percent[,5], 1), round(diff.sd.percent[,5], 1), sep = " ± ")
  df.diff$'Dry subhumid' <- paste(round(diff.prop.percent[,6], 1), round(diff.sd.percent[,6], 1), sep = " ± ")
  df.diff$'Sum drylands' <- paste(round(diff.prop.percent[,7],1), round(diff.sd.percent[,7], 1), sep = " ± ")
  df.diff$Humid <- paste(round(diff.prop.percent$Humid, 1), round(diff.sd.percent$Humid, 1), sep = " ± ")
  df.diff$Cold <- paste(round(diff.prop.percent$Cold, 1), round(diff.sd.percent$Cold, 1), sep = " ± ")
  
  k_list[[i]] <- kable(df[order(df$model, decreasing = F),], caption = i) %>% kable_styling(bootstrap_options = "bordered") %>%
    column_spec(8, italic = T, background = colorvec[5]) %>% row_spec(2, bold = T) 
  kdiff_list[[i]] <- kable(df.diff[order(df.diff$model, decreasing = F),], caption = i) %>% kable_styling(bootstrap_options = "bordered") %>% 
    column_spec(8, italic = T, background = colorvec[5]) %>% row_spec(2, bold = T) 
  
}

9.6.1 Proportion of aridity categories by continent

k_list

$AFRICA

AFRICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	27.2 ± 4.7	17.1 ± 4.9	15.6 ± 3.6	6.3 ± 1.4	66.2 ± 7.81	33.7 ± 7.1	0 ± 0
2	historical	1970_2000	26.9 ± 4.7	16.7 ± 4.2	15.1 ± 3.4	6.2 ± 1.5	64.9 ± 7.32	32.8 ± 6.2	2.3 ± 0
3	historical	1985_2015	26.5 ± 4.8	17.1 ± 4.6	15.4 ± 3.3	6 ± 1.5	65 ± 7.57	32.7 ± 6.3	2.2 ± 0
4	SSP245	2030_2060	27 ± 4.8	18 ± 4.9	16.5 ± 3.5	6.4 ± 1.7	67.9 ± 7.89	32 ± 7.6	0 ± 0
7	SSP245	2070_2100	27.9 ± 4.7	17.4 ± 5	17.2 ± 3.4	6.2 ± 1.5	68.7 ± 7.8	31.3 ± 7.9	0 ± 0
5	SSP370	2030_2060	26 ± 5	16.6 ± 3.6	17.2 ± 1.8	6.8 ± 1	66.6 ± 6.5	32.3 ± 6.7	1.1 ± 0
8	SSP370	2070_2100	28 ± 3.5	17.6 ± 4.8	16.6 ± 3.1	6.2 ± 1.6	68.4 ± 6.89	30.4 ± 6.4	1.2 ± 0
6	SSP585	2030_2060	27.1 ± 4.6	17.6 ± 4.8	16.8 ± 3.5	6.3 ± 1.5	67.8 ± 7.66	32.1 ± 7.5	0 ± 0
9	SSP585	2070_2100	29.7 ± 10.5	17.3 ± 5.3	16.2 ± 4.2	6 ± 1.5	69.2 ± 12.58	30.9 ± 8.6	0 ± 0

$ASIA

ASIA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	3.9 ± 2.3	12.2 ± 1.9	9.2 ± 1.5	4 ± 1.4	29.3 ± 3.62	32.4 ± 7.7	38 ± 9.2
2	historical	1970_2000	3.8 ± 2.1	12 ± 1.8	8.9 ± 1.5	3.8 ± 1.1	28.5 ± 3.33	30.7 ± 8.5	40.5 ± 9.6
3	historical	1985_2015	3.7 ± 2.2	11.8 ± 1.6	8.6 ± 1.5	3.7 ± 1.1	27.8 ± 3.3	30.1 ± 7.6	41.8 ± 8.8
4	SSP245	2030_2060	3.5 ± 1.8	12.5 ± 1.6	9.4 ± 1.6	4.3 ± 1.7	29.7 ± 3.35	34.3 ± 8.1	35.8 ± 9.2
7	SSP245	2070_2100	3.7 ± 2	12.5 ± 1.7	9.5 ± 1.7	4.2 ± 1.7	29.9 ± 3.56	36.2 ± 8.6	33.7 ± 9.3
5	SSP370	2030_2060	3.4 ± 2.3	12.4 ± 1.5	8.9 ± 1.7	3.7 ± 1.3	28.4 ± 3.48	34 ± 7.2	37.3 ± 9.5
8	SSP370	2070_2100	4.6 ± 1.9	12.5 ± 1.5	9.6 ± 1.5	3.8 ± 1.2	30.5 ± 3.09	35.2 ± 9.2	34.1 ± 9.7
6	SSP585	2030_2060	3.8 ± 2	12.5 ± 1.6	9.5 ± 1.8	4.2 ± 1.6	30 ± 3.52	35.4 ± 8.4	34.4 ± 9.3
9	SSP585	2070_2100	4.8 ± 3.2	13.1 ± 3.2	10.4 ± 3.4	4.2 ± 1.5	32.5 ± 5.86	36.5 ± 11.7	30.7 ± 9.9

$CENTRAL-AMERICA

CENTRAL-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	8.4 ± 3.6	30.1 ± 5.1	16.1 ± 8.7	54.6 ± 10.71	45 ± 6.1	0 ± 0
2	historical	1970_2000	0 ± 0	8.4 ± 4.1	29.7 ± 7.6	14.7 ± 7.8	52.8 ± 11.64	44.9 ± 12	2 ± 0
3	historical	1985_2015	0 ± 0	9 ± 4.2	30.7 ± 7	14.8 ± 7.4	54.5 ± 11.02	43.8 ± 11.5	1.4 ± 0
4	SSP245	2030_2060	0 ± 0	11.9 ± 6.2	36.3 ± 8	13.9 ± 6.2	62.1 ± 11.87	37.6 ± 8.6	0 ± 0
7	SSP245	2070_2100	0 ± 0	12.3 ± 5.5	38.1 ± 8.5	14.4 ± 7.9	64.8 ± 12.84	34.8 ± 8.6	0 ± 0
5	SSP370	2030_2060	0 ± 0	12.6 ± 6.3	31.4 ± 7.1	15 ± 7.9	59 ± 12.35	40.7 ± 11.8	0 ± 0
8	SSP370	2070_2100	4 ± 5.2	16 ± 5.5	35.5 ± 9.5	13.2 ± 5.2	68.7 ± 13.21	30.9 ± 10.6	0 ± 0
6	SSP585	2030_2060	0 ± 0	13 ± 7	37.6 ± 9	14.1 ± 7.1	64.7 ± 13.43	35 ± 9.4	0 ± 0
9	SSP585	2070_2100	2 ± 2.4	16.6 ± 8	41.2 ± 7.7	13.3 ± 7.7	73.1 ± 13.72	26.6 ± 8.8	0 ± 0

$EUROPE

EUROPE

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	1.9 ± 1.2	9.6 ± 4.1	4.5 ± 2.8	16 ± 5.11	49 ± 15.3	34.8 ± 20.2
2	historical	1970_2000	0 ± 0	1.5 ± 1	8.9 ± 3.4	3.9 ± 2.1	14.3 ± 4.12	44.9 ± 16.1	40.6 ± 20.4
3	historical	1985_2015	0 ± 0	1.3 ± 1	8.4 ± 2.9	3.8 ± 2	13.5 ± 3.66	46 ± 16.1	40.3 ± 19.9
4	SSP245	2030_2060	0 ± 0	1.9 ± 0.9	10.3 ± 3.6	5.1 ± 1.9	17.3 ± 4.17	52.7 ± 15.4	29.8 ± 19.7
7	SSP245	2070_2100	0 ± 0	2.2 ± 1.1	10.1 ± 4.8	5.1 ± 2.3	17.4 ± 5.44	55.4 ± 14.4	27 ± 18.6
5	SSP370	2030_2060	0 ± 0	2.1 ± 1	9.7 ± 4.7	4.8 ± 2.5	16.6 ± 5.42	52.7 ± 15.7	30.4 ± 20
8	SSP370	2070_2100	0 ± 0	3.1 ± 1.1	11.4 ± 5.9	6.4 ± 3.8	20.9 ± 7.1	52.8 ± 14.3	26 ± 18.6
6	SSP585	2030_2060	0 ± 0	2.3 ± 1.1	10.5 ± 3.7	4.9 ± 2.4	17.7 ± 4.55	55.2 ± 15.8	26.9 ± 19.9
9	SSP585	2070_2100	0 ± 0	3 ± 1.7	11.1 ± 5.4	6.1 ± 2.9	20.2 ± 6.36	57.2 ± 12.8	22.5 ± 17.3

$EUROPE-AFRICA

EUROPE-AFRICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	5.9 ± 4.9	23.6 ± 7	21.3 ± 8.8	8.8 ± 3.3	59.6 ± 12.7	24.9 ± 14.2	14.5 ± 0
2	historical	1970_2000	7 ± 5.5	26.5 ± 7.5	23.3 ± 10.3	9.9 ± 3.2	66.7 ± 14.24	27.5 ± 16.5	4.5 ± 7.5
3	historical	1985_2015	7.1 ± 5.9	26.3 ± 7.1	22.2 ± 7.7	10.3 ± 4.2	65.9 ± 12.73	27.8 ± 14.3	5.2 ± 7
4	SSP245	2030_2060	9.9 ± 6.6	28.9 ± 8.6	26.7 ± 8.9	9.4 ± 3.4	74.9 ± 14.43	23.1 ± 16.2	0.9 ± 0
7	SSP245	2070_2100	11.7 ± 5.7	28.6 ± 8.8	27.2 ± 8.1	9.5 ± 3.1	77 ± 13.61	21.9 ± 13.9	0 ± 0
5	SSP370	2030_2060	10.3 ± 6.3	28.2 ± 8.4	26 ± 8.7	9 ± 2.8	73.5 ± 13.92	24.8 ± 14	0.5 ± 0
8	SSP370	2070_2100	13.2 ± 6.2	27.3 ± 8.6	28.3 ± 9.7	8.7 ± 3.6	77.5 ± 14.81	21.2 ± 16.2	0 ± 0
6	SSP585	2030_2060	11.1 ± 6.4	28.8 ± 8.3	27.5 ± 8.8	8.7 ± 3	76.1 ± 14.01	22.6 ± 14.1	0.2 ± 0
9	SSP585	2070_2100	15.2 ± 5.7	29.1 ± 10.1	27.9 ± 8.8	8 ± 2.8	80.2 ± 14.82	18.6 ± 13.6	0 ± 0

$NORTH-AMERICA

NORTH-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	0.4 ± 0.2	5.8 ± 5.6	4.3 ± 2.6	10.5 ± 6.18	35.1 ± 9	54.2 ± 10.9
2	historical	1970_2000	0 ± 0	0.4 ± 0.2	4.3 ± 4.4	4.2 ± 2.4	8.9 ± 5.02	34.5 ± 9.8	56.5 ± 11.7
3	historical	1985_2015	0 ± 0	0.3 ± 0.2	4.6 ± 3.8	4.1 ± 2.2	9 ± 4.4	33.6 ± 10.3	57.2 ± 11.5
4	SSP245	2030_2060	0 ± 0	0.4 ± 0.3	6.5 ± 5.5	4.7 ± 2.6	11.6 ± 6.09	37.5 ± 8.7	50.8 ± 11.1
7	SSP245	2070_2100	0 ± 0	0.5 ± 0.3	6.4 ± 5.1	5.1 ± 2.9	12 ± 5.87	39.7 ± 9	48.1 ± 11
5	SSP370	2030_2060	0 ± 0	0.5 ± 0.4	5.9 ± 4.8	4.6 ± 2.7	11 ± 5.52	37 ± 9.7	51.8 ± 12.2
8	SSP370	2070_2100	0.1 ± 0	0.8 ± 0.7	7.6 ± 4	5.2 ± 2	13.7 ± 4.53	38.4 ± 10.6	47.7 ± 11.2
6	SSP585	2030_2060	0 ± 0	0.5 ± 0.4	6.6 ± 5.2	5.2 ± 2.8	12.3 ± 5.92	38.2 ± 9.1	49.3 ± 11.7
9	SSP585	2070_2100	0.2 ± 0	1 ± 1.4	7 ± 6.1	5.2 ± 2.5	13.4 ± 6.74	42.1 ± 9.3	44.5 ± 11.3

$OCEANIA

OCEANIA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	5.3 ± 0	32 ± 26	36.7 ± 17.1	6.6 ± 3.1	80.6 ± 31.27	18.9 ± 12.5	0.3 ± 0.2
2	historical	1970_2000	5.2 ± 0	31.3 ± 26.4	36.8 ± 17.4	6.8 ± 3.1	80.1 ± 31.77	19.3 ± 12.8	0.5 ± 0.4
3	historical	1985_2015	5.2 ± 0	33.4 ± 26.1	35.5 ± 17.9	6.8 ± 3.2	80.9 ± 31.81	18.4 ± 12.3	0.6 ± 0.3
4	SSP245	2030_2060	10.4 ± 0	35.8 ± 23.5	31.7 ± 15.1	5.7 ± 2.6	83.6 ± 28.05	15.9 ± 11.9	0.3 ± 0.3
7	SSP245	2070_2100	10.5 ± 0	37.9 ± 22.2	29.1 ± 13.1	6.2 ± 3.3	83.7 ± 25.99	16.2 ± 12	0.1 ± 0
5	SSP370	2030_2060	10.7 ± 0	34.8 ± 22.6	33.4 ± 15.3	5.8 ± 2.5	84.7 ± 27.41	14.9 ± 8.1	0.3 ± 0.2
8	SSP370	2070_2100	10.5 ± 10.2	38.1 ± 19.7	28.5 ± 13.8	5.7 ± 3.2	82.8 ± 26.32	16.8 ± 16.3	0.3 ± 0.3
6	SSP585	2030_2060	5.6 ± 3.6	35.6 ± 25.3	33.8 ± 16.7	6.7 ± 3.6	81.7 ± 30.74	17.7 ± 13.8	0.3 ± 0.3
9	SSP585	2070_2100	9.3 ± 6.4	39.4 ± 21.7	29.4 ± 12.9	5.5 ± 2.5	83.6 ± 26.16	16.2 ± 12.5	0 ± 0

$SOUTH-AMERICA

SOUTH-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0.1 ± 0.1	2.1 ± 1.7	15.2 ± 8.7	7.4 ± 2.6	24.8 ± 9.24	73.2 ± 11.7	1.9 ± 1
2	historical	1970_2000	0.1 ± 0.1	1.8 ± 1.3	13.9 ± 4.3	7.3 ± 2.3	23.1 ± 5.05	74.5 ± 7.3	2.3 ± 2.6
3	historical	1985_2015	0.2 ± 0.1	2 ± 1.2	14.5 ± 5.1	7.5 ± 2.4	24.2 ± 5.76	73.6 ± 7.7	2.1 ± 2.6
4	SSP245	2030_2060	0.2 ± 0.2	2.7 ± 2.2	17.4 ± 9.6	8.8 ± 3.6	29.1 ± 10.49	69.4 ± 13.1	1.3 ± 0.9
7	SSP245	2070_2100	0.2 ± 0.1	3.3 ± 3.5	18.3 ± 9.1	9.3 ± 4.6	31.1 ± 10.78	67.5 ± 14.7	1.3 ± 0.7
5	SSP370	2030_2060	0.2 ± 0.1	2.2 ± 1.7	16.3 ± 5.5	8.8 ± 3.7	27.5 ± 6.84	70.5 ± 8.9	1.9 ± 1.6
8	SSP370	2070_2100	0.2 ± 0.1	3.3 ± 2.5	18.6 ± 7	8.7 ± 3.8	30.8 ± 8.35	67.3 ± 10.8	1.8 ± 1.8
6	SSP585	2030_2060	0.2 ± 0.1	3.1 ± 2.5	18.6 ± 9	9.2 ± 4.4	31.1 ± 10.33	67.6 ± 13.1	1.2 ± 0.8
9	SSP585	2070_2100	0.6 ± 1	4.4 ± 3.8	19.6 ± 9.5	9.7 ± 3.9	34.3 ± 11	64.4 ± 12.6	1.1 ± 0.7

9.6.2 Changes in proportion of aridity categories by continent

kdiff_list

$AFRICA

AFRICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0.3 ± 6.6	0.4 ± 6.5	0.5 ± 5	0.1 ± 2.1	1.3 ± 10.7	0.9 ± 9.4	-2.3 ± 0
2	historical	1970_2000	0 ± 6.6	0 ± 5.9	0 ± 4.8	0 ± 2.1	0 ± 10.4	0 ± 8.8	0 ± 0
3	historical	1985_2015	-0.4 ± 6.7	0.4 ± 6.2	0.3 ± 4.7	-0.2 ± 2.1	0.1 ± 10.5	-0.1 ± 8.8	-0.1 ± 0
4	SSP245	2030_2060	0.1 ± 6.7	1.3 ± 6.5	1.4 ± 4.9	0.2 ± 2.3	3 ± 10.8	-0.8 ± 9.8	-2.3 ± 0
7	SSP245	2070_2100	1 ± 6.6	0.7 ± 6.5	2.1 ± 4.8	0 ± 2.1	3.8 ± 10.7	-1.5 ± 10	-2.3 ± 0
5	SSP370	2030_2060	-0.9 ± 6.9	-0.1 ± 5.5	2.1 ± 3.8	0.6 ± 1.8	1.7 ± 9.8	-0.5 ± 9.1	-1.2 ± 0
8	SSP370	2070_2100	1.1 ± 5.9	0.9 ± 6.4	1.5 ± 4.6	0 ± 2.2	3.5 ± 10.1	-2.4 ± 8.9	-1.1 ± 0
6	SSP585	2030_2060	0.2 ± 6.6	0.9 ± 6.4	1.7 ± 4.9	0.1 ± 2.1	2.9 ± 10.6	-0.7 ± 9.7	-2.3 ± 0
9	SSP585	2070_2100	2.8 ± 11.5	0.6 ± 6.8	1.1 ± 5.4	-0.2 ± 2.1	4.3 ± 14.6	-1.9 ± 10.6	-2.3 ± 0

$ASIA

ASIA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0.1 ± 3.1	0.2 ± 2.6	0.3 ± 2.1	0.2 ± 1.8	0.8 ± 4.9	1.7 ± 11.5	-2.5 ± 13.3
2	historical	1970_2000	0 ± 3	0 ± 2.5	0 ± 2.1	0 ± 1.6	0 ± 4.7	0 ± 12	0 ± 13.6
3	historical	1985_2015	-0.1 ± 3	-0.2 ± 2.4	-0.3 ± 2.1	-0.1 ± 1.6	-0.7 ± 4.7	-0.6 ± 11.4	1.3 ± 13
4	SSP245	2030_2060	-0.3 ± 2.8	0.5 ± 2.4	0.5 ± 2.2	0.5 ± 2	1.2 ± 4.7	3.6 ± 11.7	-4.7 ± 13.3
7	SSP245	2070_2100	-0.1 ± 2.9	0.5 ± 2.5	0.6 ± 2.3	0.4 ± 2	1.4 ± 4.9	5.5 ± 12.1	-6.8 ± 13.4
5	SSP370	2030_2060	-0.4 ± 3.1	0.4 ± 2.3	0 ± 2.3	-0.1 ± 1.7	-0.1 ± 4.8	3.3 ± 11.1	-3.2 ± 13.5
8	SSP370	2070_2100	0.8 ± 2.8	0.5 ± 2.3	0.7 ± 2.1	0 ± 1.6	2 ± 4.5	4.5 ± 12.5	-6.4 ± 13.6
6	SSP585	2030_2060	0 ± 2.9	0.5 ± 2.4	0.6 ± 2.3	0.4 ± 1.9	1.5 ± 4.8	4.7 ± 12	-6.1 ± 13.4
9	SSP585	2070_2100	1 ± 3.8	1.1 ± 3.7	1.5 ± 3.7	0.4 ± 1.9	4 ± 6.7	5.8 ± 14.5	-9.8 ± 13.8

$CENTRAL-AMERICA

CENTRAL-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	0 ± 5.5	0.4 ± 9.2	1.4 ± 11.7	1.8 ± 15.8	0.1 ± 13.5	-2 ± 0
2	historical	1970_2000	0 ± 0	0 ± 5.8	0 ± 10.7	0 ± 11	0 ± 16.5	0 ± 17	0 ± 0
3	historical	1985_2015	0 ± 0	0.6 ± 5.9	1 ± 10.3	0.1 ± 10.8	1.7 ± 16	-1.1 ± 16.6	-0.6 ± 0
4	SSP245	2030_2060	0 ± 0	3.5 ± 7.4	6.6 ± 11	-0.8 ± 10	9.3 ± 16.6	-7.3 ± 14.8	-2 ± 0
7	SSP245	2070_2100	0 ± 0	3.9 ± 6.9	8.4 ± 11.4	-0.3 ± 11.1	12 ± 17.3	-10.1 ± 14.8	-2 ± 0
5	SSP370	2030_2060	0 ± 0	4.2 ± 7.5	1.7 ± 10.4	0.3 ± 11.1	6.2 ± 17	-4.2 ± 16.8	-2 ± 0
8	SSP370	2070_2100	4 ± 5.2	7.6 ± 6.9	5.8 ± 12.2	-1.5 ± 9.4	15.9 ± 17.6	-14 ± 16	-2 ± 0
6	SSP585	2030_2060	0 ± 0	4.6 ± 8.1	7.9 ± 11.8	-0.6 ± 10.5	11.9 ± 17.8	-9.9 ± 15.2	-2 ± 0
9	SSP585	2070_2100	2 ± 2.4	8.2 ± 9	11.5 ± 10.8	-1.4 ± 11	20.3 ± 18	-18.3 ± 14.9	-2 ± 0

$EUROPE

EUROPE

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	0.4 ± 1.6	0.7 ± 5.3	0.6 ± 3.5	1.7 ± 6.6	4.1 ± 22.2	-5.8 ± 28.7
2	historical	1970_2000	0 ± 0	0 ± 1.4	0 ± 4.8	0 ± 3	0 ± 5.8	0 ± 22.8	0 ± 28.8
3	historical	1985_2015	0 ± 0	-0.2 ± 1.4	-0.5 ± 4.5	-0.1 ± 2.9	-0.8 ± 5.5	1.1 ± 22.8	-0.3 ± 28.5
4	SSP245	2030_2060	0 ± 0	0.4 ± 1.3	1.4 ± 5	1.2 ± 2.8	3 ± 5.9	7.8 ± 22.3	-10.8 ± 28.4
7	SSP245	2070_2100	0 ± 0	0.7 ± 1.5	1.2 ± 5.9	1.2 ± 3.1	3.1 ± 6.8	10.5 ± 21.6	-13.6 ± 27.6
5	SSP370	2030_2060	0 ± 0	0.6 ± 1.4	0.8 ± 5.8	0.9 ± 3.3	2.3 ± 6.8	7.8 ± 22.5	-10.2 ± 28.6
8	SSP370	2070_2100	0 ± 0	1.6 ± 1.5	2.5 ± 6.8	2.5 ± 4.3	6.6 ± 8.2	7.9 ± 21.5	-14.6 ± 27.6
6	SSP585	2030_2060	0 ± 0	0.8 ± 1.5	1.6 ± 5	1 ± 3.2	3.4 ± 6.1	10.3 ± 22.6	-13.7 ± 28.5
9	SSP585	2070_2100	0 ± 0	1.5 ± 2	2.2 ± 6.4	2.2 ± 3.6	5.9 ± 7.6	12.3 ± 20.6	-18.1 ± 26.7

$EUROPE-AFRICA

EUROPE-AFRICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	-1.1 ± 7.4	-2.9 ± 10.3	-2 ± 13.5	-1.1 ± 4.6	-7.1 ± 19.1	-2.6 ± 21.8	10 ± 7.5
2	historical	1970_2000	0 ± 7.8	0 ± 10.6	0 ± 14.6	0 ± 4.5	0 ± 20.1	0 ± 23.3	0 ± 10.6
3	historical	1985_2015	0.1 ± 8.1	-0.2 ± 10.3	-1.1 ± 12.9	0.4 ± 5.3	-0.8 ± 19.1	0.3 ± 21.8	0.7 ± 10.3
4	SSP245	2030_2060	2.9 ± 8.6	2.4 ± 11.4	3.4 ± 13.6	-0.5 ± 4.7	8.2 ± 20.3	-4.4 ± 23.1	-3.6 ± 7.5
7	SSP245	2070_2100	4.7 ± 7.9	2.1 ± 11.6	3.9 ± 13.1	-0.4 ± 4.5	10.3 ± 19.7	-5.6 ± 21.6	-4.5 ± 7.5
5	SSP370	2030_2060	3.3 ± 8.4	1.7 ± 11.3	2.7 ± 13.5	-0.9 ± 4.3	6.8 ± 19.9	-2.7 ± 21.6	-4 ± 7.5
8	SSP370	2070_2100	6.2 ± 8.3	0.8 ± 11.4	5 ± 14.1	-1.2 ± 4.8	10.8 ± 20.5	-6.3 ± 23.1	-4.5 ± 7.5
6	SSP585	2030_2060	4.1 ± 8.4	2.3 ± 11.2	4.2 ± 13.5	-1.2 ± 4.4	9.4 ± 20	-4.9 ± 21.7	-4.3 ± 7.5
9	SSP585	2070_2100	8.2 ± 7.9	2.6 ± 12.6	4.6 ± 13.5	-1.9 ± 4.3	13.5 ± 20.6	-8.9 ± 21.4	-4.5 ± 7.5

$NORTH-AMERICA

NORTH-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	0 ± 0.3	1.5 ± 7.1	0.1 ± 3.5	1.6 ± 8	0.6 ± 13.3	-2.3 ± 16
2	historical	1970_2000	0 ± 0	0 ± 0.3	0 ± 6.2	0 ± 3.4	0 ± 7.1	0 ± 13.9	0 ± 16.5
3	historical	1985_2015	0 ± 0	-0.1 ± 0.3	0.3 ± 5.8	-0.1 ± 3.3	0.1 ± 6.7	-0.9 ± 14.2	0.7 ± 16.4
4	SSP245	2030_2060	0 ± 0	0 ± 0.4	2.2 ± 7	0.5 ± 3.5	2.7 ± 7.9	3 ± 13.1	-5.7 ± 16.1
7	SSP245	2070_2100	0 ± 0	0.1 ± 0.4	2.1 ± 6.7	0.9 ± 3.8	3.1 ± 7.7	5.2 ± 13.3	-8.4 ± 16.1
5	SSP370	2030_2060	0 ± 0	0.1 ± 0.4	1.6 ± 6.5	0.4 ± 3.6	2.1 ± 7.5	2.5 ± 13.8	-4.7 ± 16.9
8	SSP370	2070_2100	0.1 ± 0	0.4 ± 0.7	3.3 ± 5.9	1 ± 3.1	4.8 ± 6.8	3.9 ± 14.4	-8.8 ± 16.2
6	SSP585	2030_2060	0 ± 0	0.1 ± 0.4	2.3 ± 6.8	1 ± 3.7	3.4 ± 7.8	3.7 ± 13.4	-7.2 ± 16.5
9	SSP585	2070_2100	0.2 ± 0	0.6 ± 1.4	2.7 ± 7.5	1 ± 3.5	4.5 ± 8.4	7.6 ± 13.5	-12 ± 16.3

$OCEANIA

OCEANIA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0.1 ± 0	0.7 ± 37.1	-0.1 ± 24.4	-0.2 ± 4.4	0.5 ± 44.6	-0.4 ± 17.9	-0.2 ± 0.4
2	historical	1970_2000	0 ± 0	0 ± 37.3	0 ± 24.6	0 ± 4.4	0 ± 44.9	0 ± 18.1	0 ± 0.6
3	historical	1985_2015	0 ± 0	2.1 ± 37.1	-1.3 ± 25	0 ± 4.5	0.8 ± 45	-0.9 ± 17.8	0.1 ± 0.5
4	SSP245	2030_2060	5.2 ± 0	4.5 ± 35.3	-5.1 ± 23	-1.1 ± 4	3.5 ± 42.4	-3.4 ± 17.5	-0.2 ± 0.5
7	SSP245	2070_2100	5.3 ± 0	6.6 ± 34.5	-7.7 ± 21.8	-0.6 ± 4.5	3.6 ± 41	-3.1 ± 17.5	-0.4 ± 0.4
5	SSP370	2030_2060	5.5 ± 0	3.5 ± 34.8	-3.4 ± 23.2	-1 ± 4	4.6 ± 42	-4.4 ± 15.1	-0.2 ± 0.4
8	SSP370	2070_2100	5.3 ± 10.2	6.8 ± 32.9	-8.3 ± 22.2	-1.1 ± 4.5	2.7 ± 41.3	-2.5 ± 20.7	-0.2 ± 0.5
6	SSP585	2030_2060	0.4 ± 3.6	4.3 ± 36.6	-3 ± 24.1	-0.1 ± 4.8	1.6 ± 44.2	-1.6 ± 18.8	-0.2 ± 0.5
9	SSP585	2070_2100	4.1 ± 6.4	8.1 ± 34.2	-7.4 ± 21.7	-1.3 ± 4	3.5 ± 41.2	-3.1 ± 17.9	-0.5 ± 0.4

$SOUTH-AMERICA

SOUTH-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0.1	0.3 ± 2.1	1.3 ± 9.7	0.1 ± 3.5	1.7 ± 10.5	-1.3 ± 13.8	-0.4 ± 2.8
2	historical	1970_2000	0 ± 0.1	0 ± 1.8	0 ± 6.1	0 ± 3.3	0 ± 7.1	0 ± 10.3	0 ± 3.7
3	historical	1985_2015	0.1 ± 0.1	0.2 ± 1.8	0.6 ± 6.7	0.2 ± 3.3	1.1 ± 7.7	-0.9 ± 10.6	-0.2 ± 3.7
4	SSP245	2030_2060	0.1 ± 0.2	0.9 ± 2.6	3.5 ± 10.5	1.5 ± 4.3	6 ± 11.6	-5.1 ± 15	-1 ± 2.8
7	SSP245	2070_2100	0.1 ± 0.1	1.5 ± 3.7	4.4 ± 10.1	2 ± 5.1	8 ± 11.9	-7 ± 16.4	-1 ± 2.7
5	SSP370	2030_2060	0.1 ± 0.1	0.4 ± 2.1	2.4 ± 7	1.5 ± 4.4	4.4 ± 8.5	-4 ± 11.5	-0.4 ± 3.1
8	SSP370	2070_2100	0.1 ± 0.1	1.5 ± 2.8	4.7 ± 8.2	1.4 ± 4.4	7.7 ± 9.8	-7.2 ± 13	-0.5 ± 3.2
6	SSP585	2030_2060	0.1 ± 0.1	1.3 ± 2.8	4.7 ± 10	1.9 ± 5	8 ± 11.5	-6.9 ± 15	-1.1 ± 2.7
9	SSP585	2070_2100	0.5 ± 1	2.6 ± 4	5.7 ± 10.4	2.4 ± 4.5	11.2 ± 12.1	-10.1 ± 14.6	-1.2 ± 2.7

9.6.3 Changes in proportion of aridity categories, by continent, calculated by model

tab.cont.diff <- read.table("tab.cont.diff.txt")

tab_list_percent <- list()
tab_list_sd <- list()

for(i in unique(tab.cont.diff$Continent)){
  tab.i <- subset(tab.cont.diff, Continent == i)
  
  df.percent <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "diff.mean")
  
  df.sd <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "diff.sd")
  # sd for a sum: square root of sum of squared sd
  
  drycats <- which(names(df.percent) %in% c("Hyperarid","Arid","Semi-arid","Dry subhumid"))
  
  df.percent <- df.percent %>% mutate("Sum drylands" = rowSums(.[drycats], na.rm = T))
  df.sd <- df.sd %>% mutate("Sum drylands" = round(sqrt(rowSums((.[drycats])^2, na.rm = T)),2))
  
  missing <- setdiff(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"), names(df.percent))
  
  df.percent[, missing] <- 0
  df.sd[,missing] <- 0
  
  df.percent <- df.percent %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.percent[is.na(df.percent)] <- 0
  
  tab_list_percent[[i]] <- df.percent
  
  df.sd <- df.sd %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.sd[is.na(df.sd)] <- 0
  
  tab_list_sd[[i]] <- df.sd
}

k_list <- list()

for(i in names(tab_list_percent)){
  
  df <- tab_list_percent[[i]] %>% select(c("model", "period"))
  
  df$Hyperarid <- paste(tab_list_percent[[i]]$Hyperarid, tab_list_sd[[i]]$Hyperarid, sep = " ± ")
  df$Arid <- paste(tab_list_percent[[i]]$Arid, tab_list_sd[[i]]$Arid, sep = " ± ")
  df$'Semi-Arid' <- paste(tab_list_percent[[i]][,5], tab_list_sd[[i]][,5], sep = " ± ")
  df$'Dry subhumid' <- paste(tab_list_percent[[i]][,6], tab_list_sd[[i]][,6], sep = " ± ")
  df$'Sum drylands' <- paste(tab_list_percent[[i]][,7], tab_list_sd[[i]][,7], sep = " ± ")
  df$Humid <- paste(tab_list_percent[[i]]$Humid, tab_list_sd[[i]]$Humid, sep = " ± ")
  df$Cold <- paste(tab_list_percent[[i]]$Cold, tab_list_sd[[i]]$Cold, sep = " ± ")
  
  k_list[[i]] <- kable(df[order(df$model, decreasing = F),], caption = i) %>% kable_styling(bootstrap_options = "bordered") %>%
    column_spec(8, italic = T, background = colorvec[5]) %>% row_spec(2, bold = T) 
  
}

k_list

$AFRICA

AFRICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	-5.84615384615385 ± 24.0584224816943	3 ± 22.7449628123093	4.76923076923077 ± 33.9562312547438	-0.384615384615385 ± 11.9620125225542	1.53846153846153 ± 48.91	2.84615384615385 ± 38.6972138399441	0 ± 0
2	historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
3	historical	1985_2015	-10.5384615384615 ± 19.5474439275974	11 ± 21.6024689946929	5.92307692307692 ± 11.9335768489693	-3.53846153846154 ± 10.1457329669126	2.84615384615388 ± 33.08	-2.69230769230769 ± 15.5155009339296	-2 ± 0
4	SSP245	2030_2060	-10.4615384615385 ± 41.0479706859657	24.3076923076923 ± 29.7388203514772	25.8461538461538 ± 40.2778809311971	2.53846153846154 ± 15.4844835486764	42.2307692307691 ± 66.57	-37.8461538461538 ± 53.1551912702265	0 ± 0
7	SSP245	2070_2100	10.0769230769231 ± 66.8324541153243	9 ± 46.1808040928985	42.0769230769231 ± 51.9157354733956	-2.30769230769231 ± 22.6545970882461	58.8461538461539 ± 99.03	-54.4615384615385 ± 61.6814604569523	0 ± 0
5	SSP370	2030_2060	-29.4615384615385 ± 88.4153223755319	-5.23076923076923 ± 60.9455410539522	45.4615384615385 ± 86.1293749586587	12.9230769230769 ± 45.5621581623419	23.6923076923077 ± 145	-21.3846153846154 ± 43.4291347322249	-30 ± 0
8	SSP370	2070_2100	19.4615384615385 ± 66.2578490754308	18.4615384615385 ± 58.2217247320039	31.0769230769231 ± 46.6341461779198	0.230769230769231 ± 27.395723042578	69.2307692307693 ± 103.47	-67.0769230769231 ± 52.3425600228303	-28 ± 0
6	SSP585	2030_2060	-7.84615384615385 ± 50.8524764291215	13.6153846153846 ± 39.557002038279	33.2307692307692 ± 47.566013507815	0 ± 18.4481254693623	39 ± 82.18	-34.6153846153846 ± 52.663299145829	0 ± 0
9	SSP585	2070_2100	45.5833333333333 ± 209.609792594514	9.25 ± 79.7702667323199	16.5 ± 83.830456604658	-3.83333333333333 ± 20.8232096614826	67.5 ± 240.33	-62.75 ± 98.0381419281663	0 ± 0

$ASIA

ASIA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	4.30769230769231 ± 28.0070961703893	8.76923076923077 ± 41.4993049061344	11.3076923076923 ± 61.526938023634	15.8461538461538 ± 39.580816384216	40.2307692307692 ± 88.65	90.2307692307692 ± 176.826258347073	-130.461538461538 ± 287.741091777723
2	historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
3	historical	1985_2015	-5.15384615384615 ± 19.9743425166977	-12.6153846153846 ± 41.9554097853472	-17 ± 29.8161030317511	-1.61538461538462 ± 17.1975698939243	-36.3846153846154 ± 57.83	-33.6153846153846 ± 169.076185225053	70 ± 234.12461069553
4	SSP245	2030_2060	-15.6923076923077 ± 38.0840837607712	24.2307692307692 ± 40.5301407312176	22.2307692307692 ± 88.9280175630398	27 ± 65.5222608482542	57.7692307692307 ± 123.67	192.692307692308 ± 202.871217202517	-250.461538461538 ± 375.760920308072
7	SSP245	2070_2100	-8.46153846153846 ± 31.9103310977688	23.6153846153846 ± 37.8627751702083	30.0769230769231 ± 101.227023350537	22.8461538461538 ± 72.2286256201401	68.076923076923 ± 133.85	294.461538461538 ± 208.6546650108	-362.538461538462 ± 381.791744145203
5	SSP370	2030_2060	-23.8461538461538 ± 61.3716087370566	21.6153846153846 ± 45.7156764898331	0 ± 55.9359753051052	-2.30769230769231 ± 25.9627445113975	-4.53846153846151 ± 98.28	174.230769230769 ± 120.104228239582	-169.692307692308 ± 100.362164696484
8	SSP370	2070_2100	40.0769230769231 ± 64.3343629517714	25.3076923076923 ± 63.8336178610516	32.6923076923077 ± 53.0383267323179	3.07692307692308 ± 29.0214332820347	101.153846153846 ± 108.94	241.153846153846 ± 179.945753934274	-342.307692307692 ± 156.047740886876
6	SSP585	2030_2060	-1.84615384615385 ± 30.813541811586	23.8461538461538 ± 34.0632895109632	27.0769230769231 ± 91.2044055391163	23.8461538461538 ± 54.9087821661923	72.9230769230768 ± 115.94	250 ± 201.292490338479	-322.923076923077 ± 358.702862533523
9	SSP585	2070_2100	55.0833333333333 ± 178.360441156936	59.5 ± 145.979761984636	77.9166666666667 ± 202.755768142737	22.5 ± 61.9259088089806	215 ± 313.16	324.083333333333 ± 262.604596765602	-539.083333333333 ± 416.369691136338

$CENTRAL-AMERICA

CENTRAL-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	-0.307692307692308 ± 4.67946961710789	-0.538461538461538 ± 17.9704172860815	3.30769230769231 ± 10.7964856580326	2.46153846153846 ± 21.48	-2 ± 28.2105181330179	0 ± 0
2	historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
3	historical	1985_2015	0 ± 0	1.69230769230769 ± 4.02874288466877	2.53846153846154 ± 6.25320430679856	-0.0769230769230769 ± 4.34859246311452	4.15384615384615 ± 8.62	-4 ± 5.03322295684717	-2 ± 0
4	SSP245	2030_2060	0 ± 0	9.84615384615385 ± 11.1343773500973	17.3076923076923 ± 25.7079255463648	-3.15384615384615 ± 10.0070487977738	24 ± 29.75	-23.5384615384615 ± 26.7507488512483	0 ± 0
7	SSP245	2070_2100	0 ± 0	11 ± 10.2875329080073	22.6153846153846 ± 29.2191103604543	-1.69230769230769 ± 19.0367391088242	31.9230769230769 ± 36.36	-31.4615384615385 ± 27.3484167263097	0 ± 0
5	SSP370	2030_2060	0 ± 0	11.7692307692308 ± 12.3568728929413	3.30769230769231 ± 24.015219533262	0 ± 15.6258333111123	15.0769230769231 ± 31.2	-14.6153846153846 ± 37.2615495775884	0 ± 0
8	SSP370	2070_2100	12 ± 15.556349186104	23.2307692307692 ± 15.0839530967728	18.6923076923077 ± 29.4855462200963	-3.69230769230769 ± 13.4001530798757	50.2307692307692 ± 38.97	-39.6153846153846 ± 23.8381293363471	0 ± 0
6	SSP585	2030_2060	0 ± 0	12.8461538461538 ± 14.0466256555691	21.0769230769231 ± 29.104585945808	-2.46153846153846 ± 16.4195786822489	31.4615384615384 ± 36.25	-31 ± 25.5049014897137	0 ± 0
9	SSP585	2070_2100	6 ± 7.07106781186548	25.1666666666667 ± 24.8845820588256	35.5833333333333 ± 28.5767175840541	-6.16666666666667 ± 18.5709124137396	60.5833333333333 ± 42.79	-55.0833333333333 ± 42.4273889525779	0 ± 0

$EUROPE

EUROPE

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	4.66666666666667 ± 11.211465992321	8.16666666666667 ± 13.2653570238308	7.84615384615385 ± 11.3346530151725	20.6794871794872 ± 20.74	50.3076923076923 ± 106.563271201811	-70 ± 115.036950585453
2	historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
3	historical	1985_2015	0 ± 0	-2.16666666666667 ± 9.81495457622364	-5.83333333333333 ± 24.5683955643436	-1.23076923076923 ± 8.48679215166962	-9.23076923076923 ± 27.78	12.9230769230769 ± 48.5291347859914	-4.30769230769231 ± 73.0506497979137
4	SSP245	2030_2060	0 ± 0	4.66666666666667 ± 7.41415517572279	17.5 ± 14.5695822614352	11.5 ± 9.56793888700459	33.6666666666667 ± 18.94	97.5384615384615 ± 125.997761477877	-128.538461538462 ± 146.339568233507
7	SSP245	2070_2100	0 ± 0	8.58333333333333 ± 11.1555802044076	21.6153846153846 ± 19.9354728291826	14 ± 9.32737905308882	44.1987179487179 ± 24.68	123.076923076923 ± 120.805257569405	-166.615384615385 ± 139.523318517932
5	SSP370	2030_2060	0 ± 0	7.16666666666667 ± 8.48349556313298	17.6153846153846 ± 11.2215452110249	10.3846153846154 ± 7.83728759987005	35.1666666666667 ± 16.1	90.5384615384615 ± 75.058438771195	-125.153846153846 ± 73.5026599902413
8	SSP370	2070_2100	0 ± 0	19.25 ± 13.7452471124319	37.6923076923077 ± 28.0338623078847	30 ± 31.8093277724213	86.9423076923077 ± 44.57	92.2307692307692 ± 110.037988778235	-177.692307692308 ± 82.173986775403
6	SSP585	2030_2060	0 ± 0	9.83333333333333 ± 8.20014781832768	19.0833333333333 ± 14.2028059540957	12.2307692307692 ± 8.41777728138339	41.1474358974358 ± 18.43	125.923076923077 ± 134.35987839782	-164.846153846154 ± 144.399241776545
9	SSP585	2070_2100	0 ± 0	18.2727272727273 ± 18.9530520449394	37.4166666666667 ± 24.7403179715487	26.0833333333333 ± 24.6298353197659	81.7727272727273 ± 39.72	150.666666666667 ± 131.040821768245	-230.916666666667 ± 148.445556798893

$EUROPE-AFRICA

EUROPE-AFRICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	-2.07692307692308 ± 5.1552810861216	-1.15384615384615 ± 9.04476048923864	1.76923076923077 ± 14.8893353677156	-0.384615384615385 ± 8.18065259762795	-1.84615384615385 ± 19.92	1 ± 20.7163381577601	13 ± 0
2	historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
3	historical	1985_2015	0.538461538461538 ± 5.82545255202524	0.230769230769231 ± 8.23765587767384	-4.15384615384615 ± 17.8505189179762	2 ± 13.1782649338472	-1.38461538461538 ± 24.37	2.46153846153846 ± 16.7163362444017	-6.5 ± 9.19238815542512
4	SSP245	2030_2060	11.2307692307692 ± 8.13586551586908	7 ± 10.0249688278817	11.6153846153846 ± 19.8432963556061	-3.30769230769231 ± 10.1438373359117	26.5384615384615 ± 25.76	-22.2307692307692 ± 35.5109228410871	-54 ± 0
7	SSP245	2070_2100	19.3846153846154 ± 8.20100056434642	6.84615384615385 ± 11.3860013016434	14.9230769230769 ± 23.0704050624082	-2.53846153846154 ± 10.6584503612187	38.6153846153846 ± 29.03	-36.5 ± 35.8570393748183	0 ± 0
5	SSP370	2030_2060	14 ± 6.7700320038633	6.07692307692308 ± 9.92019437361267	10.5384615384615 ± 14.3097599829358	-4.53846153846154 ± 8.04793332286189	26.076923076923 ± 20.34	-23.0833333333333 ± 20.7121322958894	-56 ± 0
8	SSP370	2070_2100	25.3076923076923 ± 16.7101995568805	-0.923076923076923 ± 23.063179668256	17.6923076923077 ± 22.1260955110499	-6.53846153846154 ± 21.33313301188	35.5384615384615 ± 41.9	-30.9230769230769 ± 28.5115810927815	0 ± 0
6	SSP585	2030_2060	16.9230769230769 ± 9.8696634395635	7.84615384615385 ± 10.1722347089693	16.5384615384615 ± 17.1736978381447	-5.84615384615385 ± 13.0054229320321	35.4615384615384 ± 25.79	-33.1666666666667 ± 32.0137281159003	-57 ± 0
9	SSP585	2070_2100	32.4166666666667 ± 12.273092026005	11.3333333333333 ± 19.90583894847	21.1666666666667 ± 22.0735684528448	-8.66666666666667 ± 18.1174618598308	56.25 ± 36.91	-55.5454545454545 ± 27.4494576863137	0 ± 0

$NORTH-AMERICA

NORTH-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	0.375 ± 3.06768875307034	31.75 ± 66.2901679258422	4.92307692307692 ± 27.7172074665467	37.0480769230769 ± 71.92	21.6923076923077 ± 91.341287319759	-56.0769230769231 ± 173.993036995958
2	historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
3	historical	1985_2015	0 ± 0	0.333333333333333 ± 4.18330013267038	-3 ± 32.2038958682161	-2.69230769230769 ± 36.415163085782	-5.35897435897436 ± 48.79	-19.7692307692308 ± 70.0180379689808	25.0769230769231 ± 113.976943237409
4	SSP245	2030_2060	0 ± 0	3.36363636363636 ± 2.87307247638229	60 ± 104.176132903207	14.1538461538462 ± 34.693530025655	77.5174825174826 ± 109.84	82.1538461538462 ± 94.51705856074	-159.153846153846 ± 214.335735919859
7	SSP245	2070_2100	0 ± 0	4.7 ± 3.05686840482943	56.9230769230769 ± 85.9946331062405	25.4615384615385 ± 49.6548342470556	87.0846153846154 ± 99.35	145.923076923077 ± 122.363162715515	-231.923076923077 ± 215.561847868333
5	SSP370	2030_2060	0 ± 0	5.45454545454545 ± 6.89004552036685	44.4615384615385 ± 45.1582686865787	12.0769230769231 ± 33.2452040111591	61.9930069930071 ± 56.5	70.4615384615385 ± 48.5397009066039	-131.615384615385 ± 77.0427786595846
8	SSP370	2070_2100	4 ± 0	15.6923076923077 ± 20.0307455984736	91.5384615384615 ± 59.6093049009065	28.3846153846154 ± 49.8339550600098	139.615384615385 ± 80.24	108 ± 110.324521299664	-243.923076923077 ± 105.99564577414
6	SSP585	2030_2060	0 ± 0	7.09090909090909 ± 7.86707754448303	62.0769230769231 ± 90.9537075829068	29.8461538461538 ± 46.234630155772	99.013986013986 ± 102.33	103.384615384615 ± 97.9187915754159	-201.307692307692 ± 203.199157074771
9	SSP585	2070_2100	5 ± 0	20.1818181818182 ± 38.4468937154048	91.1666666666667 ± 135.959107042812	30.75 ± 35.4481054757258	147.098484848485 ± 145.67	198.75 ± 99.7835156726801	-339.583333333333 ± 220.157158054023

$OCEANIA

OCEANIA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	1 ± 0	6 ± 42.7668095606862	-1.23076923076923 ± 38.0660694191775	-1.23076923076923 ± 6.80874249468831	4.53846153846154 ± 57.66	-3 ± 12.0623380818148	0.166666666666667 ± 1.16904519445001
2	historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
3	historical	1985_2015	1 ± 0	5.08333333333333 ± 32.8977985815944	-2.23076923076923 ± 21.2922280271223	1.76923076923077 ± 9.42650382480028	5.62179487179487 ± 40.3	-3 ± 7.85281265959316	-0.25 ± 0.957427107756338
4	SSP245	2030_2060	45 ± 0	50.1538461538462 ± 62.0991762611042	-26.4615384615385 ± 39.9408216085903	-5.84615384615385 ± 12.4354744839522	62.8461538461539 ± 74.87	-19.8461538461538 ± 18.3296033501644	-1.5 ± 0.707106781186548
7	SSP245	2070_2100	45 ± 0	66.6153846153846 ± 81.5041292016734	-48.4615384615385 ± 58.5172558376521	-2.46153846153846 ± 13.0103509087661	60.6923076923076 ± 101.18	-17.3846153846154 ± 18.8172724800845	-5 ± 0
5	SSP370	2030_2060	47 ± 0	42.3846153846154 ± 50.2535877020047	-12.3846153846154 ± 33.3105050035232	-5.30769230769231 ± 10.3068311925038	71.6923076923077 ± 61.17	-27 ± 61.1882341631134	-3.33333333333333 ± 3.21455025366432
8	SSP370	2070_2100	70.6666666666667 ± 84.0376899571456	63.1538461538462 ± 110.812789690425	-57 ± 62.8887907341205	-6.84615384615385 ± 13.452671072109	69.974358974359 ± 153.22	-14.1538461538462 ± 57.5135435554765	-5 ± 0
6	SSP585	2030_2060	23.5 ± 0.707106781186548	33.8461538461538 ± 55.8910340958878	-23.2307692307692 ± 43.1164196220609	-0.769230769230769 ± 12.2348807796524	33.3461538461538 ± 71.65	-12 ± 20.0956048262632	-1.5 ± 0.707106781186548
9	SSP585	2070_2100	60.3333333333333 ± 48.6757982300582	69.8333333333333 ± 126.826892282594	-50 ± 91.7645207731576	-9 ± 14.1806526327567	71.1666666666666 ± 164.55	-23.9166666666667 ± 22.0802599347897	0 ± 0

$SOUTH-AMERICA

SOUTH-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	-0.142857142857143 ± 0.377964473009227	5.07692307692308 ± 13.901208211648	20.0769230769231 ± 83.3551253557747	1.30769230769231 ± 18.8121619145019	26.3186813186813 ± 86.58	-20.6153846153846 ± 79.5115698725019	-6.81818181818182 ± 35.0679859182651
2	historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
3	historical	1985_2015	0.5 ± 1.22474487139159	2.69230769230769 ± 5.87912430248779	10.5384615384615 ± 15.6504280272425	3 ± 9.31844049899624	16.7307692307692 ± 19.18	-14.2307692307692 ± 13.1539211273917	-2.54545454545455 ± 5.75088925929958
4	SSP245	2030_2060	1.57142857142857 ± 1.13389341902768	13.1538461538462 ± 26.6328310982459	55.6153846153846 ± 93.1374776531432	23.4615384615385 ± 24.7810928754679	93.8021978021979 ± 100	-78.4615384615385 ± 97.2167812885335	-18.8 ± 42.023273974734
7	SSP245	2070_2100	1.42857142857143 ± 0.975900072948533	22.6153846153846 ± 49.907478500285	68.9230769230769 ± 90.7647706422684	31 ± 41.3481156362254	123.967032967033 ± 111.53	-107 ± 130.065368180773	-22.7777777777778 ± 43.6886077192264
5	SSP370	2030_2060	1.25 ± 0.886405260427918	5.92307692307692 ± 19.7250329043306	38.1538461538462 ± 31.8821950149979	24.0769230769231 ± 36.5637833984704	69.4038461538462 ± 52.38	-59.4615384615385 ± 55.635143846037	-12.8888888888889 ± 18.6711304186734
8	SSP370	2070_2100	2 ± 2.26778683805536	22.9230769230769 ± 23.8168761541809	73.5384615384615 ± 82.0869613932032	22.6923076923077 ± 35.8849657827727	121.153846153846 ± 92.73	-109.461538461538 ± 109.343964461248	-15 ± 15.1492574075431
6	SSP585	2030_2060	1.5 ± 1.0690449676497	19.2307692307692 ± 29.9197644992789	72.8461538461538 ± 84.535836733942	30.1538461538462 ± 38.093844983685	123.730769230769 ± 97.44	-107.538461538462 ± 99.3576162024628	-20.1 ± 41.1675438503036
9	SSP585	2070_2100	7.375 ± 15.259072242908	40.5 ± 51.561966250681	96.9166666666667 ± 122.811354474571	38 ± 40.7274756488451	182.791666666667 ± 140.12	-160.166666666667 ± 133.095204020781	-29.375 ± 47.7850469737732

9.6.4 Changes in proportion of aridity categories, figure

list_plots <- list()

for(i in names(tab_list_percent)){
  df.percent <- tab_list_percent[[i]]
  df.sd <- tab_list_sd[[i]]
  
  df <- merge(df.percent, df.sd, by = c("model", "period"))
  
  g <- ggplot(data = df)+geom_point(aes(x=period,y= get("Sum drylands.x"), col = model), shape = "\u2605", size = 5, position = position_dodge(width = 0.5))+
    geom_errorbar(aes(x = period, ymin = get("Sum drylands.x")-get("Sum drylands.y"), ymax = get("Sum drylands.x")+get("Sum drylands.y"), col = model), position = position_dodge(width = 0.5))+
    scale_color_manual(values = colorvec[c(1,3,9,11)])+
    ylim(0,100)+
    labs(x="", y = "Percent of drylands", title = i)+
    theme_minimal()
  list_plots[[i]] <- g
}

ggarrange(plotlist = list_plots, ncol = 2, nrow = 4, common.legend = T, legend = "bottom")

tab.flow <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","POLAR","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, Continent, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model, Continent) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL)

df245 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP245")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+
  labs(title = "SSP 2-4.5", y = "%", x = "")+
  theme_minimal()

df370 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP370")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+  labs(title = "SSP 3-7.0", y = "%", x = "")+
  theme_minimal()


df585 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP585")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+  labs(title = "SSP 5-8.5", y = "%", x = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

10 Changes of variables over time and scenarios

10.1 Temperature

quantile(cmip6s$t.mean, probs = seq(0,1,0.1)) 
##          0%         10%         20%         30%         40%         50% 
## -53.4587313 -42.5828261 -28.4644864 -17.7502317  -7.2080554  -0.2580962 
##         60%         70%         80%         90%        100% 
##   7.1198406  15.3318920  22.8072968  26.4971129  34.3727997

t.breaks <- c(-50,-40,-30,-20,-10,0,10,20,30)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")

map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = t.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = t.breaks, palette = function(x) colscale,
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Mean annual temperature, °C")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = t.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = t.breaks, palette = function(x) colscale,
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Mean annual temperature, °C", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

#labels = c("2030-2060", "2070-2100"

10.2 Precipitation

quantile(cmip6s$pr.mean, probs = seq(0,1,0.1)) 
##          0%         10%         20%         30%         40%         50% 
##    8.579696   55.000319  104.132310  200.852004  337.042710  463.946055 
##         60%         70%         80%         90%        100% 
##  579.741473  712.099280  959.490392 1431.315602 6863.723388

p.breaks <- c(10,50,100,200,300,400,500,600,700)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")

map_list <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = pr.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = p.breaks, palette = function(x) rev(colscale),
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Annual mean precipitation")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = pr.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = p.breaks, palette = function(x) rev(colscale),
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Annual mean\nprecipitation, mm", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")